Dekompozycja klasyczna - estymacja trendu i sezonowości - Grupa Żywiec S.A.

Time Series Przychodów Grupy Żywiec

# Importuje Time Series zawierający dane przychodow Grupy Zywiec
load("/Users/wspolne/daneTS.RData")

# importowane dane to Time Series przychodów    GŻ podzielone na kwartały - zgodnie z czestotliwością ich raportowania przez spółkę

dane.ts

##        Qtr1   Qtr2   Qtr3   Qtr4
## 2005  61057 149585  53639  76316
## 2006  61379 194508  62531  89736
## 2007  94247 276643  -5537  74607
## 2008 105011 242350   1379  80031
## 2009  79331 181047  63547  74982
## 2010  81931 185390  60206  82205
## 2011  85264 190220  72367  81266
## 2012  80006 188721  65808  90017
## 2013  86143 174774  68797  88478
## 2014  86785 165124  64649  80355
## 2015  86808 168711  74626  93995
## 2016  97119 186404  80502  96460
## 2017 104067 194324  89953 107316
## 2018 112729 214280  99393 126014
## 2019 125194 240354 106869 118456
## 2020 140697 246861 132527 153084
## 2021 160422 305486 130196 179544

# tak natomiast wygląda wykres linowy zaimportowanych danych
plot(dane.ts)

Na powyższym wykresie, zauważalne jest występowanie wielkości determinsitycznych(nielosowych), tj. sezowoności, jak tendencji długoterminowej (trendu)./ Aby lepiej to zobrazować, posłuże się dwoma wykresami z pakietu forecast

Wykres kwartalny i sezonowy

# importuje bibliotekę forecast
library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

# wykres kwartalny

monthplot(dane.ts, main="wykres kwartalny")

# wykres sezonowy
seasonplot(dane.ts, year.labels = TRUE, col=rainbow(3),main="wykres sezonowy")

Autokorelacja

ACF

Istotną kwestią, wiążącą się również z analizą szeregów czasowych jest występowanie autokorelacji, czyli korelacji danych w czasie. Do analizy występowania tego zjawiska posłuży wykres funkcji autokorelacji ACF.

Wykres ten pokazuje stopień podobieństwa między danym szeregiem czasowym a jego opóźnioną wersją, przesuniętą o parametr opóźnienia h jednostek czasowych.

# wykres autokorelacji ACF
acf(dane.ts, main="wykres ACF")

Na powyższym wykresie, widoczne jest występowanie sezonowości, poprzez oscylacje okresowe o okresie 4, odpowiadjące sezonowości kwartalnej.

PACF

Pojęciem, które również związane jest z funckją ACF jest również funkcja autokorelacji cząsteczkowej PACF. PACF(h) mierzy zależność pomiędzy obserwacjami odległymi o h jednostek czasowych

pacf(dane.ts, main="wykres PACF")

Widoczna na powyższym wykresie wysoka, dodatnia wartość opóźnienia dla h=1 może świadczyć o trendu wzrostowego.

Dekompozycja

Potwierdzenie występowania w danych zarówno sezonowości, jak i trendu długoterminowego uzyskać można również poprzez dekompozycje szeregu czasowego, przy użyciu funkcji decompose(). Przy jej użyciu, wyznaczana jest dekompozycja danego szeregu na elementy składowe tj. trend, sezonowość oraz losowe wahania. Z uwagi na charakter danych, dekompozycja będzie podjęta dla modelu addytywnego

# dekompozycja addytywna
dane.ts.dekomp.add <- decompose(dane.ts, "additive")
plot(dane.ts.dekomp.add)

Uzyskany model, potwierdza występowanie zarówno trendu długoterminowego (o wyraźnym zintensyfikowaniu po roku 2015) oraz sezonowości (cykliczne wahania o powtarzalnym interwale występowania).

Dekompozycja na podstawie modelu regresji

Metoda dekompozycji na podstawie modelu regresji, jest metodą parametryczną, tak więc dzięki jej zastosowaniu jesteśmy w stanie utrzymać funkcyjną postać trendu m(t) Funkcja m(t) może przyjmować postać:
liniową m(t)=a+bt
kwadratową m(t)=a+bt+ct^2
wykładniczą m(t)=a*exp(bt)
lub inną.

Funkcja tslm()

funkcja tslm() z pakietu forecast służy do odpowiedniego dopasowania modeli liniowych dla szeregów czasowych zapisannych w klasie Time Series. Funkcja, umożliwia również wykorzystanie zmiennych trend oraz season w celu zdefiniowania występowania zjawisk trendu oraz sezonowości

library(forecast)
# wykorzystuje funkcję tslm z wraz ze zmiennymi trend oraz season
dane.tslm.trend.sezonowosc <- tslm(dane.ts ~trend + season)
# podstawowe własności dopasowanego modelu
summary(dane.tslm.trend.sezonowosc)

## 
## Call:
## tslm(formula = dane.ts ~ trend + season)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -50325 -11893  -3449   9588  97541 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  59742.0     8126.1   7.352 4.86e-10 ***
## trend         1127.6      158.6   7.109 1.29e-09 ***
## season2     108083.7     8792.2  12.293  < 2e-16 ***
## season3     -27357.4     8796.5  -3.110  0.00281 ** 
## season4       -755.0     8803.6  -0.086  0.93193    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25630 on 63 degrees of freedom
## Multiple R-squared:  0.8385, Adjusted R-squared:  0.8282 
## F-statistic: 81.76 on 4 and 63 DF,  p-value: < 2.2e-16

jakość dopasowania powyższego modelu obrazuje wartość Adjusted R-squared: wynosząca: 0.8282. Miara ta przedstawia skorygowany współczynnik determinacji - Identyfikuje ona jaki procent zmienności w zmiennej zależnej(objaśnianej), jest wyjaśniony przez zbudowany model.

# zestawienie wykresu przychodów Grupy Żywiec oraz zbudowanego modelu
plot(dane.ts, main="dekompozycja na podstawie regresji liniowej", col="black")
lines(fitted(dane.tslm.trend.sezonowosc), col="red",lty=3)
legend("topleft",legend = c("przychody Grupy Zywiec","dopasowany model"), col=c("black","red"), lty=c(1,3))

Tak natomiast wygląda graficzne zestawienie danych oraz zbudowanego modelu.

Analiza reszt z modelu dekompozycji

Następnym krokiem podjętym w ocenie jakości zbudowanego modelu jest analiza reszt z modelu dekompozycji

# wykres reszt modelu opartego na regresji liniowej
tsdisplay(residuals(dane.tslm.trend.sezonowosc))

Powyższy wykres pokazuje, że dopasowanie reszt nie jest zadowalające - na początku analizowanego modelu, wariancja reszt jest wysoka, oraz dodatkowo - widoczne są pozostałości efektów sezonowych. Problem ten może wynikać ze specyfiki danych tj. wzrostu wahań sezonowych wraz z poziomem szeregu.
Jednym z potencjalnych rozwiązań tego problemu jest użycie transformacji Boxa-Coxa - w tym transformacji logarytmicznej, która powinna umożliwić skuteczniejsze dopasowanie modelu addytywnego.

Transformacje Boxa-Coxa - transformacja logarytmiczna

Pośród transformacji Boxa-Coxa, jednym z przekszłaceń wstępnych, które znajduję swoje zastosownie przy modelowaniu danych sezonowych jest transformacja logarytmiczna.
Logartymując dane wejściowe, możliwe jest ustablizowanie wariancji oraz zniwelowanie dużego wpływu wartości odstających. Transformację te można zastosowować poprzez zdefiniowanie odpowiedniej wartości zmiennej lambda w funkcji tslm(). Jako, że zastosowana zostanie transformacja logarytmiczna, to wartość lambda=0

library(forecast)
# dekompozycja danych zlogarytmowanych, transformacja logarytmiczna z lambda=0

dane.log.tslm <- tslm(dane.ts ~trend + season, lambda=0)

## Warning in log(x): NaNs produced

# podstawowe własności dopasowanego nowego modelu
summary(dane.log.tslm)

## 
## Call:
## tslm(formula = dane.ts ~ trend + season, lambda = 0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5194 -0.0968  0.0175  0.1635  0.6145 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.03453    0.15721  70.188  < 2e-16 ***
## trend        0.01259    0.00309   4.076 0.000133 ***
## season2      0.75549    0.16926   4.463 3.47e-05 ***
## season3     -0.47490    0.17220  -2.758 0.007634 ** 
## season4     -0.01303    0.16949  -0.077 0.938960    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4934 on 62 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5208, Adjusted R-squared:  0.4899 
## F-statistic: 16.85 on 4 and 62 DF,  p-value: 2.134e-09

# wykres przychodów Grupy Żywiec oraz utworzonego modelu za pomocą transformacji logarytmicznej

plot(dane.ts, main="dekompozycja na podstawie regresji liniowej - transformacja logarytmiczna", col="black")
lines(fitted(dane.log.tslm), col="blue",lty=3)
legend("topleft",legend = c("przychody Grupy Zywiec","nowy dopasowany model"), col=c("black","blue"), lty=c(1,3))

Warto podkreślić, że z uwagi na zmianę skali, współczynnik determinacji R^2 nie może być porównywany ze współczynnikami, otrzymanymi w oryginalnym modelu, przed zastosowaniem przekształcenia.

Analiza reszt z modelu dekompozycji po przekształceniu

W celu analizy jakości otrzymanego modelu, należy przyjrzeć się resztom z modelu dekompozycji.

# wykres reszt modelu dekompozycji po zastosowaniu przekształcenia 
tsdisplay(residuals(dane.log.tslm))

Na powyższym wykres widoczny pozytywny jest efekt zastosowania transformacji logarytmicznej, jej użycie pozwoliło na usunięcie z modelu multiplikatywnego charakteru sezonowości oraz znaczego wpływu wartości odstających.
Aby być pewnym co do losowości reszt nowgo modelu, należy przeprowadzić test Ljung-box.

Test Ljung-Boxa

Test Ljung-Boxa jest testem statystycznym stosowanym w analizie szeregów czasowych do sprawdzenia, czy szereg czasowy charakteryzuje się brakiem autokorelacji, będąc tzw. “białym szumem”.

Autokorelacja oznacza, że wartości w szeregu czasowym są skorelowane z wartościami z przeszłości. Oznacza to, że dane mają tendencję do powtarzania okresowych wzorców. Biały szum natomiast oznacza brak autokorelacji w danych. Oznacza to, że każda wartość w szeregu czasowym jest losowa.

Najczęściej test Ljunga-Box jest stosowany do badania autokorelacji w resztach modelu, a nie w samym szeregu czasowym.

H0=reszty modelu są losowe
H1=reszty modelu nie są losowe

Przyjęty poziom istotności(p-value) to 0.05

# test Ljung-Boxa dla reszt modelu po zastosowaniu przekształcenia logarytmicznego
Box.test(residuals(dane.log.tslm),type="Lj",lag=4)

## 
##  Box-Ljung test
## 
## data:  residuals(dane.log.tslm)
## X-squared = 0.81375, df = 4, p-value = 0.9366

Otrzymana wartość p=0.9366, tak więc jest to więcej, niż przyjęty poziom istotności. Zatem nie mamy podstaw do odrzucenia H0 i możemy stwierdzić, że reszty modelu są losowe.

Prognozy oparte na dekompozycji

TBA