[1] 56488 54867 50274 58463 67610 69213Szeregi czasowe – projekt zaliczeniowy
1. Temat projektu
Dane przedstawiają liczbę fabrycznie nowych samochodów osobowych i ciągników siodłowych zarejestrowanych po raz pierwszy na terytorium Polski. Są to dane kwartalne pochodzące z lat 2010 - 2022.
[1] 1035 1324 1527 2396 2917 3129 Qtr1 Qtr2 Qtr3 Qtr4
2010 56488 54867 50274 58463
2011 67610 69213 63622 74148 Qtr1 Qtr2 Qtr3 Qtr4
2010 1035 1324 1527 2396
2011 2917 3129 2562 3160
2. Identyfikacja nielosowych składowych
2.1. Trend
Do aproksymacji trendu w szeregach za pomocą wielomianu użyjemy poniższych funkcji:
2.1.1. Samochody osobowe
Wybór stopnia wielomianu dla samochodów osobowych:
Wybieramy wielomian stopnia 4.
2.1.2. Ciągniki siodłowe
Wybór stopnia wielomianu dla ciągników siodłowych:
Wybieramy wielomian stopnia 5.
2.1.3. Zbudowanie modeli
2.2. Sezonowość
2.2.1 Samochody osobowe
Na powyższych korelogramach przedstawiono funkcję autokorelacji (ACF). Korelacje nie występują po różnicowaniu szeregu, natomiast są widoczne przed różnicowaniem. Wskazuje to na obecność sezonowości. Wniosek ten jest zgodny z charakterem danych, które prezentują sezonowość związaną z kwartalnymi okresami czasowymi.
Augmented Dickey-Fuller Test
data: diff(X)
Dickey-Fuller = -4.948, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary
KPSS Test for Level Stationarity
data: diff(X)
KPSS Level = 0.15348, Truncation lag parameter = 3, p-value = 0.1Po różnicowaniu szereg liczby samochodów osobowych jest stacjonarny, ponieważ wyeliminowaliśmy lokalne trendy.
2.2.2 Ciągniki siodłowe
Na podstawie powyższego korelogramu dla drugiego szeregu możemy wywnioskować, że tak jak w przypadku szeregu pierwszego, obecna jest sezonowość.
Augmented Dickey-Fuller Test
data: diff(Y)
Dickey-Fuller = -3.7075, Lag order = 3, p-value = 0.03271
alternative hypothesis: stationary
KPSS Test for Level Stationarity
data: diff(Y)
KPSS Level = 0.075118, Truncation lag parameter = 3, p-value = 0.1Szereg liczby ciągników siodłowych jest, po różnicowaniu, również stacjonarny.
3. Identyfikacja (modelowanie) reszt
3.1. Korelacje
Do zbadania korelacji reszt w modelach użyjemy testu Breuscha-Godfreya.
\(H_0\): Błędy nie są skorelowane
\(H_1\): Błędy są skorelowane
Durbin-Watson test
data: model
DW = 1.2675, p-value = 0.0003344
alternative hypothesis: true autocorrelation is greater than 0Dla pierwszego modelu wartość \(p\) jest mniejsza niż przyjęty poziom istotności \(0.05\), zatem odrzucamy hipotezę zerową.
Durbin-Watson test
data: model2
DW = 1.6223, p-value = 0.01478
alternative hypothesis: true autocorrelation is greater than 0Dla drugiego modelu również odrzucamy hipotezę o braku korelacji reszt.
3.2. Stacjonarność
Do zbadania stacjonarności szeregów użyjemy rozszerzonego testu Dickey’a-Fullera i testu Kwiatkowskiego-Phillipsa-Schmidta-Shina. Dla testu ADF formułujemy hipotezy:
\(H_0\): Szereg czasowy jest niestacjonarny.
\(H_1\): Szereg czasowy jest stacjonarny.
Natomiast dla testu KPSS:
\(H_0\): Szereg czasowy jest stacjonarny
\(H_1\): Szereg czasowy jest niestacjonarny
Augmented Dickey-Fuller Test
data: X
Dickey-Fuller = -1.0175, Lag order = 3, p-value = 0.9271
alternative hypothesis: stationary
KPSS Test for Level Stationarity
data: X
KPSS Level = 1.0592, Truncation lag parameter = 3, p-value = 0.01
Augmented Dickey-Fuller Test
data: Y
Dickey-Fuller = -2.5026, Lag order = 3, p-value = 0.3722
alternative hypothesis: stationary
KPSS Test for Level Stationarity
data: Y
KPSS Level = 1.0508, Truncation lag parameter = 3, p-value = 0.01Na podstawie powyższych testów możemy stwierdzić, że oba szeregi nie są stacjonarne.
3.3. Homoskedastyczność
Użyjemy testu Goldfielda-Quandta do zbadania homoskedastyczności.
\(H_0\): Wariancja każdego błędu w modelu jest taka sama.
\(H_1\): Błędy mają różne wariancje
Goldfeld-Quandt test
data: model
GQ = 3.7084, df1 = 21, df2 = 21, p-value = 0.002042
alternative hypothesis: variance increases from segment 1 to 2
Goldfeld-Quandt test
data: model2
GQ = 2.1804, df1 = 20, df2 = 20, p-value = 0.04451
alternative hypothesis: variance increases from segment 1 to 2Zarówno dla pierwszego jak i drugiego modelu odrzucamy hipotezę zerową na korzyść hipotezy alternatywnej. Oznacza to, że wariancje są niejednorodne.
3.4. GARCH
Modele mają niejednorodną wariancję i nieskorelowane reszty, zatem wybieramy dopasowanie modelem GARCH.
Dla obu szeregów czasowych zastosowaliśmy model z parametrami (1,3), ponieważ przy przeszukiwaniu różnych parametrów, te okazały się mieć najmniejszy współczynnik AIC.
Call:
garch(x = X, order = c(1, 3))
Model:
GARCH(1,3)
Residuals:
Min 1Q Median 3Q Max
0.4392 0.7375 0.8124 0.8662 1.0299
Coefficient(s):
Estimate Std. Error t value Pr(>|t|)
a0 6.437e+08 NA NA NA
a1 5.548e-01 NA NA NA
a2 5.180e-01 NA NA NA
a3 5.189e-01 NA NA NA
b1 1.266e-06 NA NA NA
Diagnostic Tests:
Jarque Bera Test
data: Residuals
X-squared = 9.0511, df = 2, p-value = 0.01083
Box-Ljung test
data: Squared.Residuals
X-squared = 1.4861, df = 1, p-value = 0.2228Wartość \(p\) w teście Jarque-Bera jest mniejsza niż \(0.05\), zatem odrzucamy hipotezę o normalności rozkładu reszt.
Na podstawie testu Boxa-Ljunga nie mamy podstaw do odrzucenia hipotezy o braku autokorelacji.
Call:
garch(x = Y, order = c(1, 3))
Model:
GARCH(1,3)
Residuals:
Min 1Q Median 3Q Max
0.2425 0.6972 0.8110 0.9367 1.3091
Coefficient(s):
Estimate Std. Error t value Pr(>|t|)
a0 2.178e+06 NA NA NA
a1 5.634e-01 NA NA NA
a2 5.343e-01 NA NA NA
a3 4.633e-01 NA NA NA
b1 2.208e-08 NA NA NA
Diagnostic Tests:
Jarque Bera Test
data: Residuals
X-squared = 1.7329, df = 2, p-value = 0.4204
Box-Ljung test
data: Squared.Residuals
X-squared = 0.73164, df = 1, p-value = 0.3924Dla drugiego modelu wartość \(p\) w teście Jarque-Bera jest większa niż przyjęty poziom istotności, zatem nie mamy podstaw do odrzucenia hipotezy o normalności rozkładu.
Na podstawie testu Boxa-Ljunga stwierdzamy, tak jak dla pierwszego modelu, że reszty nie są skorelowane.
4. SARIMA
Series: ts_data
ARIMA(0,1,2)
Coefficients:
ma1 ma2
-0.2171 -0.2799
s.e. 0.1327 0.1176
sigma^2 = 181186520: log likelihood = -556.36
AIC=1118.72 AICc=1119.23 BIC=1124.51Do zbudowania modelu SARIMA użyjemy parametrów dla składowych niesezonowych \((0,1,2)\) oraz parametrów sezonowych \((P,D,Q)\) gdzie:
\(P=4\) - stopień integracji sezonowej
\(D=1\) - stopień różnicowania sezonowego
\(Q=2\) - stopień ruchomej średniej sezonowej
Series: ts_data
ARIMA(0,1,2)(4,1,2)[4]
Coefficients:
ma1 ma2 sar1 sar2 sar3 sar4 sma1 sma2
-0.0718 -0.4801 0.7830 -0.0029 -0.1928 -0.1097 -1.9648 0.9846
s.e. 0.1454 0.1693 0.4681 0.1839 0.2276 0.2397 1.2936 1.2970
sigma^2 = 137584358: log likelihood = -514.48
AIC=1046.97 AICc=1051.83 BIC=1063.62
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -1414.245 10158.16 6298.331 -2.13377 6.996225 0.4335062
ACF1
Training set -0.09615871
Series: ts_data2
ARIMA(1,1,0)
Coefficients:
ar1
-0.4911
s.e. 0.1250
sigma^2 = 1e+06: log likelihood = -424.29
AIC=852.59 AICc=852.84 BIC=856.45W drugim przypadku używamy parametrów dla składowych niesezonowych \((1,1,0)\). Parametry sezonowe pozostają bez zmian.
Series: ts_data2
ARIMA(1,1,0)(4,1,2)[4]
Coefficients:
ar1 sar1 sar2 sar3 sar4 sma1 sma2
-0.2697 0.7805 0.0613 -0.1734 -0.0563 -1.9214 0.9935
s.e. 0.1463 0.2455 0.2259 0.2055 0.3097 1.0069 1.0272
sigma^2 = 846209: log likelihood = -390.58
AIC=797.16 AICc=800.95 BIC=811.96
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -59.04279 806.8021 538.5454 -6.456574 16.20475 0.5109741
ACF1
Training set -0.03740173
Możemy zauważyć, że lepiej dopasowany do danych jest model drugi (dotyczący ciągników siodłowych), ponieważ wartości Kryterium Informacyjnego Akaikego (AIC) i Bayesowskiego Kryterium Schwarza (BIC) są niższe niż w przypadku modelu pierwszego (dotyczącego samochodów osobowych).
5. Holt-Winters
Skorzystaliśmy z funkcji `HoltWinters` , w której wskazaliśmy że nasz model jest multiplikatywny, ponieważ wartości przez cały czas rosną, tym samym rośnie rozmiar sezonowych wahań.
5.1 Samochody osobowe
Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
Call:
HoltWinters(x = ts_data, seasonal = "multiplicative")
Smoothing parameters:
alpha: 0.7137117
beta : 0.03224566
gamma: 0.3049768
Coefficients:
[,1]
a 1.079146e+05
b 1.041016e+03
s1 9.892742e-01
s2 9.840274e-01
s3 9.219543e-01
s4 9.765842e-01Parametr \(alpha\) informuje nas o tym, jak duży wpływ na zmianę szeregu czasowego mają bieżące obserwacje. Jest on równy \(0.71\), co świadczy o tym, że ten wpływ jest duży.
Na podstawie parametru \(beta\) równego \(0.03\) możemy stwierdzić, że ostatnie obserwacje mają bardzo mały wpływ na zmianę trendu.
Parametr \(gamma\) wynosi \(0.3\), co oznacza że ostatnie obserwacje mają mały wpływ na zmiany wzorców sezonowych.
Przez ostatnie spadki i stabilizacje w danych, przedział ufności dla przyszłych danych jest szerszy, zakładając zarówno mocny spadek jak i wzrost w sprzedaży aut osobowych.
Dla porównania poniżej przedstawiamy ten sam model, jednak dla danych bez ostatnich 5 lat (do 2018 roku):
Widać wyraźnie, że model przewiduje wzrost danych, zachowując przy tym sezonowe spadki.
Teraz pokażemy wykres reszt dla modelu Holt-Winters i korelogram reszt.
Linie przerywane reprezentują przedziały ufności dla autokorelacji. Wartości przekraczające te przedziały sugerują istotne korelacje.
5.2 Ciągniki siodłowe
Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
Call:
HoltWinters(x = ts_data2, seasonal = "multiplicative")
Smoothing parameters:
alpha: 0.7260936
beta : 0.04199913
gamma: 0.1959948
Coefficients:
[,1]
a 7805.5256524
b 176.3652696
s1 0.9362797
s2 1.0217929
s3 0.8560945
s4 1.0227962Tak jak w przypadku samochodów osobowych, bieżące obserwacje mają duży wpływ na zmianę szeregu czasowego \((alpha=0.73)\), ostatnie obserwacje mają bardzo mały wpływ na zmianę trendu \((beta=0.04)\) oraz mają one mały wpływ na zmiany wzorców sezonowych \((gamma=0.2)\).
Przedziały ufności mają duże zakresy. Wyniki prognoz wskazują na tendencję wzrostową liczby ciągników siodłowych.
