Cvičenie 9: SARIMA modelovanie
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 9
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.09837353Zopakujte 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$coefCvič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$coef3 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.250004 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.79Keďž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.62Keďž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.3176004Keďž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.003059Sezó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.63Vytvorte 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:
