1 Inštalácia knižníc v R-ku

Začneme s inštaláciou knižníc potrebných na toto cvičenie. Spustením nasledovného príkazu sa nainštalujú tie z knižníc, ktoré nainštalované nie sú, a zároveň sa aj načítajú.

# názvy balíkov, ktoré budú potrebné pre cvičenie 5
packages <- c("astsa", "urca", "datasets", "ggplot2", "seasonal","fpp3", "climatol")

# inštalácia balíkov, ak nie sú nainštalované
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages],dependencies=TRUE)
}

# načítanie balíkov
invisible(lapply(packages, library, character.only = TRUE))

2 Zápis modelov v R-ku

Uvažujme najprv klasický ARMA(1,1) proces, napríklad \((1 - 0.7 L)x_t = (1 + 0.8 L)u_t\), ktorý možno prepísať na \(x_t = 0.7 x_{t-1} + u_t + 0.8 u_{t-1}\).

V R-ku vieme takýto proces simulovať pomocou funkcie arima.sim ako

set.seed(123)
x <- arima.sim(model = list(ar = 0.7, ma = 0.8), n = 1000)

Pre simulovaný proces x vieme potom spätne odhadnúť ARMA koeficienty.

m <- sarima(x, 1, 0, 1, details = FALSE)
m$fit$coef
##        ar1        ma1      xmean 
## 0.67321585 0.79735653 0.09837353

Zopakujte uvedený postup pre sezónne procesy.

Cvičenie 1: Sezónny AR proces \((1 - 0.8 L)(1 + 0.9 L^4) x_t = u_t\).

x <- arima.sim(...)
m <- sarima(...)
m$fit$coef

Cvičenie 2: Sezónny MA proces \(x_t = (1 - 0.8 L)(1 + 0.9 L^4) u_t\)

x <- arima.sim(...)
m <- sarima(...)
m$fit$coef

3 Stacionarita a invertovateľnosť sezónnych modelov

Overovanie stacionarity (AR časť) a invertovateľnosti (MA časť) sa robí analýzou koreňov príslušných charakteristických polynómov. Pri sezónnych modeloch môžeme skúmať nesezónne a sezónne polynómy buď osobitne, alebo rozpísať celý model do jedného polynómu. Oba postupy vedú k rovnakým koreňom. Proces je stacionárny, ak všetky korene AR polynómu ležia mimo jednotkovej kružnice \((|z| > 1)\). Proces je invertovateľný, ak všetky korene MA polynómu ležia mimo jednotkovej kružnice \((|z| > 1)\).

Cvičenie: Overte stacionaritu sezónneho AR procesu \((1 - 0.8 L)(1 + 0.9 L^4) x_t = u_t\).

Mali by ste dostať nasledujúce absolútne hodnoty koreňov:

## [1] 1.02669 1.02669 1.02669 1.02669 1.25000

4 Modelovanie sezónnych časových radov

Hlavnou náplňou tohto cvičenia je modelovanie sezónnych časových radov. Budeme hľadať SARIMA model pre rôzne reálne dáta.

4.1 Mesačné teploty z balíka climatol

Budeme pracovať s dátami z balíka climatol, ktorý sa využíva na analýzu a vizualizáciu klimatologických a meteorologických dát. Konktrétne budeme pracovať a dátmi Temp.dat. Tieto dáta obsahujú mesačné teploty v piatich staniciach počas obdobia 1961-2005. V tomto cvičení budeme pracovať s dátami z tretej stanice, t. j. Temp.dat[,3].

data(climatol_data)
data <- Temp.dat[,3]
data_ts <- ts(data, start=c(1961,1), frequency = 12)

Na modelovanie použijeme dáta x z rokov 1961 - 2000, teda bez posledných 5 rokov. Dáta z rokov 2001-2005 uložíme do premennej x.pred a použijeme ich na zhodnotenie kvality predikcií.

x <- window(data_ts,  end = c(2000,12), frequency = 12)
x.pred <-  window(data_ts, start=c(2001,1))

Dáta vykazujú jasnú sezónnosť, ich stredná hodnota je niekde okolo 17, dáta nevykazujú zjavný lineárny trend. Testujeme preto prítomnosť jednotkového koreňa pomocou funkcie ur.df s prametrami lags=15 a type="drift".

summary(ur.df(x, type="drift", lags = 15, selectlags = "BIC" ))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3921 -0.6956 -0.0300  0.7027  2.7056 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.91982    1.60119   2.448 0.014744 *  
## z.lag.1      -0.21851    0.08947  -2.442 0.014978 *  
## z.diff.lag1  -0.32384    0.09508  -3.406 0.000719 ***
## z.diff.lag2  -0.35917    0.09260  -3.879 0.000121 ***
## z.diff.lag3  -0.44403    0.08746  -5.077 5.63e-07 ***
## z.diff.lag4  -0.56496    0.07884  -7.166 3.19e-12 ***
## z.diff.lag5  -0.57055    0.07238  -7.882 2.44e-14 ***
## z.diff.lag6  -0.46855    0.06855  -6.835 2.68e-11 ***
## z.diff.lag7  -0.51413    0.06382  -8.055 7.17e-15 ***
## z.diff.lag8  -0.64082    0.05706 -11.231  < 2e-16 ***
## z.diff.lag9  -0.68810    0.05265 -13.069  < 2e-16 ***
## z.diff.lag10 -0.48943    0.05323  -9.194  < 2e-16 ***
## z.diff.lag11 -0.32579    0.05107  -6.380 4.41e-10 ***
## z.diff.lag12 -0.13544    0.04674  -2.898 0.003942 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.053 on 450 degrees of freedom
## Multiple R-squared:  0.8711, Adjusted R-squared:  0.8674 
## F-statistic:   234 on 13 and 450 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -2.4423 3.008 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.44 -2.87 -2.57
## phi1  6.47  4.61  3.79

Keďže testová štatistika (-2.4423) je väčšia ako kritická hodnota (-2.87), hypotézu o jednotkovom koreni nezamietame. Naviac, ako sme predpokladali, v regresii sa odhadli lagy 1-12, preto diferencujeme sezónne.

Priebeh sezónne diferencovaných dát:

V dátach po sezónnom diferencovaní nevidíme ani trend ani sezónnosť. Naviac, keďže pôvodné dáta nemali trend, pri teste jednotkového koreňa pre diferencie volíme type="none". Dostávame:

summary(ur.df(x.d12, type="none", lags = 10, selectlags = "BIC" ))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7937 -0.8722 -0.0427  0.8687  4.7920 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.75896    0.05729 -13.248   <2e-16 ***
## z.diff.lag  0.01837    0.04680   0.393    0.695    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.329 on 455 degrees of freedom
## Multiple R-squared:  0.3731, Adjusted R-squared:  0.3703 
## F-statistic: 135.4 on 2 and 455 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -13.2484 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Keďže tentokrát je testová štatistika menšia ako kritická hodnota, hypotézu o jednotkovom koreni zamietame. Sezónne diferencované dáta sú stacionárne a ďalej ich netreba diferencovať.

Zobrazíme si výberovú ACF a PACF:

acf2(x.d12)

Pozrime sa najprv na hodnoty ACF a PACF pre celé roky. Hodnota ACF pre lag 12 je výrazne nenulová, pre lag 24 je na hranici, pre lagy 36,48 sú autokorelácie takmer nulové. PACF je nenulová pre všetky lagy 12, 24, 36, 48. To nás motivuje skúsiť pri modelovaní použiť jeden až dva sezónne MA členy.

Aby sme zachytili autokorelácie a parciálne autokorelácie na začiatku, musíme do modelu pridať klasické členy. Priebeh však nie je jednoznačný - zdá sa, že ani ACF ani PACF sa po niektorom lagu nevynuluje. Obe majú naviac výrazne nenulové prvé hodnoty. Skúsme preto použiť zmiešaný ARMA(1,1) model. Teda spolu s jednym sezónnym MA členom, uvažujme model sarima(x, 1, 0, 1, 0, 1, 1, 12). Dostávame

sarima(x, 1, 0, 1, 0, 1, 1, 12)   
## initial  value 0.315762 
## iter   2 value 0.099550
## iter   3 value 0.060789
## iter   4 value 0.044382
## iter   5 value 0.032343
## iter   6 value 0.030805
## iter   7 value 0.029043
## iter   8 value 0.028678
## iter   9 value 0.028349
## iter  10 value 0.028312
## iter  11 value 0.028298
## iter  12 value 0.028055
## iter  13 value 0.028038
## iter  14 value 0.028036
## iter  15 value 0.028035
## iter  16 value 0.028035
## iter  17 value 0.028035
## iter  18 value 0.028035
## iter  18 value 0.028035
## iter  18 value 0.028035
## final  value 0.028035 
## converged
## initial  value 0.013209 
## iter   2 value 0.009007
## iter   3 value 0.001575
## iter   4 value 0.000763
## iter   5 value 0.000627
## iter   6 value 0.000453
## iter   7 value 0.000230
## iter   8 value -0.000066
## iter   9 value -0.000244
## iter  10 value -0.000425
## iter  11 value -0.000443
## iter  12 value -0.000457
## iter  13 value -0.000460
## iter  14 value -0.000460
## iter  15 value -0.000460
## iter  16 value -0.000461
## iter  17 value -0.000461
## iter  18 value -0.000461
## iter  19 value -0.000461
## iter  19 value -0.000461
## final  value -0.000461 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##          Estimate     SE  t.value p.value
## ar1        0.4179 0.1327   3.1482  0.0017
## ma1       -0.1118 0.1455  -0.7687  0.4425
## sma1      -0.9660 0.0634 -15.2250  0.0000
## constant   0.0025 0.0005   4.6737  0.0000
## 
## sigma^2 estimated as 0.9334676 on 464 degrees of freedom 
##  
## AIC = 2.858323  AICc = 2.858507  BIC = 2.902644 
##

Všetky p-hodnoty Ljung-Boxovho testu sú väčšie ako 5%. Model vyzerá dobre. Z tabuľky koeficientov však vidíme, že koeficient ma1 nie je signifikantný (p-value=0.4425), preto ho skúsime vynechať.

model.100.011 <- sarima(x, 1, 0, 0, 0, 1, 1, 12)   
## initial  value 0.315762 
## iter   2 value 0.067943
## iter   3 value 0.044568
## iter   4 value 0.035001
## iter   5 value 0.032564
## iter   6 value 0.028697
## iter   7 value 0.028119
## iter   8 value 0.028089
## iter   9 value 0.028068
## iter  10 value 0.028059
## iter  11 value 0.028059
## iter  12 value 0.028059
## iter  13 value 0.028059
## iter  13 value 0.028059
## iter  13 value 0.028059
## final  value 0.028059 
## converged
## initial  value 0.013003 
## iter   2 value 0.008473
## iter   3 value 0.001225
## iter   4 value 0.000496
## iter   5 value 0.000168
## iter   6 value 0.000154
## iter   7 value 0.000144
## iter   8 value 0.000140
## iter   9 value 0.000140
## iter  10 value 0.000140
## iter  11 value 0.000140
## iter  12 value 0.000140
## iter  13 value 0.000140
## iter  13 value 0.000140
## iter  13 value 0.000140
## final  value 0.000140 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##          Estimate     SE  t.value p.value
## ar1        0.3176 0.0440   7.2171       0
## sma1      -0.9640 0.0622 -15.5036       0
## constant   0.0025 0.0005   4.8005       0
## 
## sigma^2 estimated as 0.9356825 on 465 degrees of freedom 
##  
## AIC = 2.855251  AICc = 2.855362  BIC = 2.890708 
##

Všetky p-hodnoty Ljung-Boxovho testu sú stále väčšie ako 5%. Model vyzerá dobre. Hodnota BIC sa v porovnaní s predošlým modelom znížila z hodnoty 2.902644 na hodnotu 2.890708. Všetky koeficienty sú signifikantné.

Overme stacionaritu a invertovateľnosť. Či je proces stacionárny závisí od AR časti:

koef <- model.100.011$fit$coef
koef_ar <- koef[grep("^ar", names(koef))]
koef_ar
##       ar1 
## 0.3176004

Keďže jediný AR koeficient je menší ako 1, proces je stacionárny.

Či je proces invertovateľný závisí od (sezónnej) MA časti:

koef_sma <- koef[grep("^sma", names(koef))]
abs(polyroot(c(1,rep(0,11), koef_sma)))
##  [1] 1.003059 1.003059 1.003059 1.003059 1.003059 1.003059 1.003059 1.003059
##  [9] 1.003059 1.003059 1.003059 1.003059

Sezónne MA korene sú mimo jednotkového kruhu, teda model je invertovateľný.

Na záver urobme predikcie z tohto modelu na 5 rokov a porovnajme ich so skutočnými hodnotami, ktoré sme si uložili do premennej x.pred.

predikcie_sarima<-sarima.for(x, 60, 1, 0, 0, 0, 1, 1, 12,details = FALSE)
points(x.pred, col = "blue", type = "l")

4.2 Dáta exp o čínskom exporte

Pozrime sa ďalej na dáta exp o čínskom exporte z knižnice seasonal.

plot(exp)

Aby sme odstránili rastúcosť disperzie, budeme pracovať s logaritmami dát. Na modelovanie použijeme dáta x z rokov 2000 - 2012, teda bez posledného roku. Dáta za posledný rok uložíme do premennej x.pred a použijeme ich na zhodnotenie kvality predikcií.

x <- window(..., start= ... , end = ... )
x.pred <-  window(..., start=...)
plot(x)

Zodpovedajte otázky:

  • Sú dáta stacionárne?
  • Je v dátach prítomný trend?
  • Pozorujeme v dátach sezónnosť?
  • Je potrebné dáta diferencovať? Ak áno, sezónne alebo klasicky? Koľkokrát?

Pre stacionárne dáta, ktoré ďalej netreba diferencovať, by ste mali dostať nasledovnú ACF a PACF.

acf2(x.d12.d1, main = "")

  • Aké tipy na sarima modely máte na základe ACF a PACF?

  • Určte najlepší z vyhovujúcich modelov na základe BIC.

  • Pre najlepší z modelov spravte predikcie na jeden rok. Porovnajte s reálnymi dátami odloženými v premennej x.pred.

Výsledok pre porovnanie:

4.3 Harmonizovaný index spotrebiteľských cien na Slovensku

Analyzujme harmonizovaný index spotrebiteľských cien (HICP) na Slovensku. HICP je štatistický ukazovateľ používaný na meranie inflácie, resp. zmien v cenovej hladine spotrebiteľských tovarov a služieb. Tento index je štandardizovaný pre krajiny Európskej únie a ďalšie štáty, aby umožnil porovnanie medzi krajinami.

Zdroj dát: https://fred.stlouisfed.org/series/CP0000SKM086NEST

Dáta môžete do R-ka jednoducho načítať príkazom

HICP.SK <- read.table("https://anna-hlubinova.github.io/HICP_Slovakia.txt", header = TRUE)
head(HICP.SK)
##       DATE Value
## 1 1.1.1996 43.01
## 2 1.2.1996 43.12
## 3 1.3.1996 43.22
## 4 1.4.1996 43.35
## 5 1.5.1996 43.56
## 6 1.6.1996 43.63

Vytvorte z dát časový rad. Dáta za rok 2024 odložte do premennej hicp.pred. Pomocou týchto dát zhodnotíme predikcie. Pri modelovaní budeme pracovať so skráteným časovým radom do konca roka 2023.

hicp.ts <- ts(... , frequency = ..., start = ...)

hicp <- window(hicp.ts, end = ...)
hicp.pred <- window(hicp.ts , start = ...)

Priebeh časového radu:

plot(hicp, type = "l")

Zodpovedajte otázky:

  • Sú dáta stacionárne?
  • Je v dátach prítomný trend?
  • Pozorujeme v dátach sezónnosť?
  • Je potrebné dáta diferencovať? Ak áno, sezónne alebo klasicky? Koľkokrát?

Pre stacionárne dáta, ktoré ďalej netreba diferencovať, by ste mali dostať nasledovnú ACF a PACF.

  • Aké tipy na sarima model máte na základe ACF a PACF?

  • Určte najlepší z vyhovujúcich modelov na základe BIC.

  • Pre najlepší z modelov spravte predikcie na jeden rok. Porovnajte s reálnymi dátami odloženými v premennej hicp.pred.

Výsledok pre porovnanie:

4.4 Dáta o nezamestnanosti v USA

Modelujme dáta o nezamestanosti v USA dostupné v balíku fpp3 pod názvom us_employment. Informácie o dátach nájdete v helpe pomocou ?us_employment.

Konkrétne budeme modelovať dáta v sektore Leisure and Hospitality (rekreácia a pohostinstvo) od januára 2001 do septembra 2019. Z dát vytvoríme časový rad.

leisure <- us_employment %>%
  filter(Title == "Leisure and Hospitality",
         year(Month) > 2000) %>%
  mutate(Employed = Employed/1000) %>%
  select(Month, Employed)

leisure <- ts(leisure$Employed, start = c(2001,1), end = c(2019,9),
              frequency = 12)

Do premennej leisure.pred si odložíme posledný rok dát na zhodnotenie predikcií. Pri modelovaní budeme preto pracovať s časovým radom skráteným o tento jeden rok.

leisure.pred <- window(leisure, start = ... )
leisure <- window(leisure, end = ... )

Priebeh časového radu:

Zodpovedajte otázky:

  • Sú dáta stacionárne?
  • Je v dátach prítomný trend?
  • Pozorujeme v dátach sezónnosť?
  • Je potrebné dáta diferencovať? Ak áno, sezónne alebo klasicky? Koľkokrát?

Pre stacionárne dáta, ktoré už netreba diferencovať, by ste mali dostať nasledovnú ACF a PACF.

  • Aké tipy na sarima modely máte na základe ACF a PACF?

  • Určte najlepší z vyhovujúcich modelov na základe BIC.

  • Pre najlepší z modelov spravte predikcie na jeden rok. Porovnajte s reálnymi dátami odloženými v premennej leisure.pred.

Výsledok pre porovnanie:

5 Modelovanie trendu

5.1 Exponenciálne zhladzovanie

  • Majme dáta \(x_1, x_2, \dots x_n\) a chceme predikovať hodnotu \(x_{n+k}\)
  • Predpoklad: v dátach nie je trend ani sezónnosť
  • Model: \[ x_t = \mu_t + w_t ,\] kde \(\mu_t\) je stredná hodnota (môže závisieť od času) a \(w_t\) sú nezávislé náhodné odchýlky s nulovou strednou hodnotu.
  • Označme \(a_t\) náš odhad strednej hodnoty \(\mu_t\)
  • Základná myšlienka: ďalší odhad strednej hodnoty bude váženým priemerom predchádzajúceho odhadu (teda \(a_{t-1}\)) a novej relizovanej hodnoty \(x_t\) : \[ a_t = \alpha x_t + (1-\alpha) a_{t-1}, \] kde \(a \in (0,1)\) je parameter zhladzovania.
  • Platí:
    • \(a \approx 1\): slabé zhladzovanie, \(a_t \approx x_t\),
    • \(a \approx 0\): silné zhladzovanie, \(a_t \approx a_{t-1}\).
  • Keďže nemáme trend ani sezónnosť, vieme spraviť iba \[ \hat{x}_{n+k|n} = a_n \]
  • Pre daný parameter \(\alpha\):
    • \(a_1 = x_1\) a potom rekurentne
    • máme teda predikčné chyby: \(e_t = x_t - \hat{x}_{t|t-1} = a_n\)
  • Optimálny parameter: minimalizuje sumu štvorcov predikčných chýb: \[ \sum_{t=2}^n e_t \rightarrow \min_{\alpha} \]

5.1.1 Dáta o prietokoch Nílu

Postup aplikujeme na dáta o prietokoch Nílu z Cvičenia 2.

Pozn: Dáta neobsahujú sezónnosť, pretože ide o ročné údaje.

x <- Nile
plot(x)

V R-ku použijeme všeobecnejšiu funkciu HoltWinters, pričom zadáme parametre beta=FALSE, gamma=FALSE, čo zodpovedá práve exponenciálnemu zhladzovaniu.

model <- HoltWinters(x, beta=FALSE, gamma=FALSE)
model
## Holt-Winters exponential smoothing without trend and without seasonal component.
## 
## Call:
## HoltWinters(x = x, beta = FALSE, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.2465579
##  beta : FALSE
##  gamma: FALSE
## 
## Coefficients:
##       [,1]
## a 805.0389

Vykreslenie:

plot(model)

Hodnotu sum of squared errors, podľa ktorej sa vybrala optimálna \(\alpha\) vieme vypísať pomocou príkazu:

model$SSE
## [1] 2038872

Otázka: Viete nájsť fixné \(\alpha\) také, že dostanete menšiu hodnotu SSE?

Cvičenie: Vykreslite závislosť SSE od parametra \(\alpha\) pre tieto dáta. Výpočet má potvrdiť optimálnu hodnotu, ktorú našlo R-ko.

Výstup na porovnanie:

Detail:

5.2 Holt-Wintersova metóda

Charakteristiky časového radu:

  • \(a_t\) = level, sezónne očistená stredná hodnota

  • \(b_t\) = slope, zmena hodnoty level z jednej periódy na druhú (zachytáva rôzne, aj krátkodobé trendy)

  • \(s_t\) = seasonal component, sezónna zložka (závisí napr. od mesiaca)

Typ sezónnosti:

  • aditívna - napr. v januári roku \(i\) je hodnota o 100 vyššia ako v januári roku \(i-1\)
  • multiplikatívna - napr. v januári roku \(i\) je hodnota o 10 percent vyššia ako v januári roku \(i-1\)

5.2.1 Dáta exp o čínskom exporte

Vráťme sa k dátam exp o čínskom exporte z knižnice seasonal. Na modelovanie opäť použijeme dáta z rokov 2000-2012, teda vynecháme posledný rok. Ten použijeme na zhodnotenie kvality predikcií.

x <- window(log(exp), start=c(2000,1), end = c(2012,12))
x.pred <-  window(log(exp), start=c(2013,1))

plot(x)

Cvičenie 1: Odhadnite optimálne parametre Holt-Wintersovej metódy (s aditívnou sezónnosťou, keďže sme logaritmovaním dosiahli stabilizovanie disperzie) a spravte predikcie na rok 2013.

Na nájdenie predikcií použite funkciu predict. Predikcie spolu s reálnymi dátami môžete vykresliť napríklad pomocou funkcie ts.plot.

Predikcie:

Cvičenie 2: Spravte predikcie z Holt-Wintersovej metódy s aditívnou sezónnosťou a pevne zvolenými parametrami \(\alpha=\beta=\gamma=\frac{1}{2}\) (teda každý člen sa berie s rovnakou váhou).

Predikcie:

Vykreslite do jedného obrázka dáta, predikcie nájdené pomocou SARIMA modelu a obe predikcie nájdené pomocou Holt-Wintersovej metódy.