1. Wybór szeregu czasowego i jego opis

Wybrany zbiór danych reprezentuje średnią, miesięczną ilość pasażerów lotnicznych pewngo lotu, na przełomie lat 1946 - 1960.

Analizie poddane zostają 144 kolejne obserwacje, w odstępie miesięcznym. Jest to szereg czasowy okresów.

Źródło: https://www.kaggle.com/datasets/ashfakyeafi/air-passenger-data-for-time-series-analysis/

2. Analiza wstępna wybranego szeregu czasowego

- Ładowanie zbioru danych

- Tworzenie obiektu typu ts

- Podgląd surowych danych liczbowych

- Wykres szeregu czasowego

data <- read.csv("AirPassengers.csv", header = TRUE)
ts <- ts(data$X.Passengers, start = c(1949, 1), frequency = 12)

print(ts)
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432
plot(ts, main = "Miesięczna liczba pasażerów lotnicztch w latach 1949 - 1960", xlab = "Czas", ylab = "Liczba pasażerów")

- Analiza anomalii i brakujących wartości

anyNA(ts)
## [1] FALSE
boxplot(ts, main = "Boxplot liczby pasażerów", col = "lightblue")

Nie znaleziono żadnych brakujących wartości lub “outlierów”.

- Analiza zmienności i wskaźnik zmienności (CV)

mean_value <- mean(ts)
std_dev <- sd(ts)
cat("Średnia liczba pasażerów:", mean_value, "\n")
## Średnia liczba pasażerów: 280.2986
cat("Odchylenie standardowe liczby pasażerów:", std_dev, "\n")
## Odchylenie standardowe liczby pasażerów: 119.9663
cv <- (std_dev / mean_value) * 100
cat("Współczynnik zmienności (CV):", cv, "%\n")
## Współczynnik zmienności (CV): 42.79947 %

- Analiza dynamiki

ts_diff <- diff(ts)
plot(ts_diff, type = "o", col = "purple", main = "Dynamika liczby pasażerów", 
     xlab = "Czas", ylab = "Zmiana liczby pasażerów")

Odchylenie standardowe informuje nas o rozproszeniu wartości liczby pasażerów względem średniej. W tym przypadku wynosi ono około 120. Oznacza to, że liczba pasażerów w poszczególnych miesiącach może różnić się od średniej o około 120 pasażerów, co pokazuje znaczną zmienność w danych.

Współczynnik zmienności (CV) jest stosunkiem odchylenia standardowego do średniej, wyrażony w procentach. Wynik 42.8% oznacza, że zmienność liczby pasażerów w analizowanym okresie jest stosunkowo wysoka, bo odchylenie standardowe stanowi niemal 43% średniej wartości. Może to świadczyć o dużych wahaniach w liczbie pasażerów pomiędzy poszczególnymi miesiącami, co może być efektem sezonowości, czynników ekonomicznych, czy innych zmiennych.

Wykres zmiany liczby pasażerów (dynamiki) pokazuje różnice między kolejnymi obserwacjami szeregu czasowego. Rysując te różnice, możemy dostrzec zmiany w dynamice liczby pasażerów, np. w okresach wzrostów lub spadków. Zmienność na wykresie może również sugerować, czy w danych występują wyraźne zmiany sezonowe, nagłe skoki lub spadki liczby pasażerów, które są istotne do dalszej analizy.

Podsumowanie:

  1. Wysoka zmienność: Współczynnik zmienności na poziomie 42.8% wskazuje na to, że liczba pasażerów jest dość zmienna w analizowanym okresie.

  2. Sezonowość i inne czynniki: Duża zmienność może sugerować wpływ czynników takich jak sezonowość, zmieniające się warunki ekonomiczne, zmiany cen biletów, a także wydarzenia zewnętrzne (np. wojny, zmiany polityczne). Warto rozważyć dalszą dekompozycję szeregu czasowego, aby oddzielić komponenty sezonowe, cykliczne i losowe.

- Określenie stacjonarności

adf.test(ts)
## Warning in adf.test(ts): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts
## Dickey-Fuller = -7.3186, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(ts)
## Warning in kpss.test(ts): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts
## KPSS Level = 2.7395, Truncation lag parameter = 4, p-value = 0.01

Test ADF wskazuje, że szereg jest stacjonarny (p-value = 0.01, hipoteza alternatywna = szereg stacjonarny)

Test KPSS wskazuje, że szereg NIE jest stacjonarny (p-value = 0.01, hipoteza alternatywna = szereg niestacjonarny)

W związku z tym, testy zostają ponowione na zróżniocowanym zbiorze.

adf.test(ts_diff)
## Warning in adf.test(ts_diff): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_diff
## Dickey-Fuller = -7.0177, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(ts_diff)
## Warning in kpss.test(ts_diff): p-value greater than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_diff
## KPSS Level = 0.014626, Truncation lag parameter = 4, p-value = 0.1

Przy zróżnicowanym szeregu, oba testy wskazują, że szereg jest stacjonarny.

ADF: p-value = 0.01, hipoteza alternatywna = szereg stacjonarny,

KPSS: p_value = 0.1, hipoteza alternatywna = szereg niestacjonarny

W związku z tym, przy predykcji przyszłych wartości szeregu wykorzystany zostanie szereg zróżnicowany.

3. Dekompozycja Szeregu Czasowego

decomp <- decompose(ts)
plot(decomp)

Porównując poszczególne składowe szeregu wnioskować można, że szereg jest addytywny. Dostrzec można liniową (proporcjonalną) zależność wartości poszczególnej składowej od wartości całego szeregu, co wskazuje na charakter addytywny.

Automatyczna dekopozycja skutecznie wydzieliła 3 składowe szeregu: trend, seasonal, random

4. Wybrna Analiza zaawansowana

- Analiza wykresów ACF i PACF celem ustalenia parametrów do modelu ARIMA

acf(ts_diff, main = "Autokorelacja (ACF)", lag.max = 30)

pacf(ts_diff, main = "Częściowa autokorelacja (PACF)", lag.max = 30)

Na podstaiwe powyższych wykresów można ustalić parametry:

  • AR (p): Wykres PACF pokazuje znaczący skok przy lag 1, co sugeruje komponent autoregresyjny rzędu 1 (p=1)

  • I (d): Wynosi 0, ponieważ dane są już stacjonarne.

  • MA (q): Wykres ACF pokazuje znaczące skoki przy lag 1 oraz 2, co sugeruje komponent średniej ruchomej rzędu 2 (q=2)

  • Seasonal AR (p): brak istotnych skoków w wartościach sezonowych na wykresie PACF co sugeruje wartość 0

  • Seasonal I (d): dane mają wyraźny wzorzec sezonowy, ale nie są stacjonarne w sensie sezonowym, należy zastosować sezonowe różnicowanie (d=1)

  • Seasonal MA (q): brak istotnych skoków w wartościach sezonowych na wykresie ACF co sugeruje wartość 0

Sezonowy model ARIMA(1,0,2)(0,1,0)[12]

model <- Arima(ts_diff, order = c(1, 0, 2), seasonal = c(0, 1, 0))
summary(model)
## Series: ts_diff 
## ARIMA(1,0,2)(0,1,0)[12] 
## 
## Coefficients:
##           ar1     ma1      ma2
##       -0.4862  0.1661  -0.0854
## s.e.   0.4709  0.4698   0.1956
## 
## sigma^2 = 140:  log likelihood = -508.08
## AIC=1024.16   AICc=1024.48   BIC=1035.66
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE    MASE        ACF1
## Training set 0.2400279 11.19328 8.272793 NaN  Inf 0.87823 0.004177268
forecasted_diff <- forecast(model, h = 12)
plot(forecasted_diff)

Predykcja została wykonana na szeregu zróżnicowanym (stacjonarnym).

Wnioski:

  1. Analiza szeregu czasowego pozwoliła zidentyfikować ogólny trend wzrostowy liczby pasażerów.

  2. Model predykcyjny wykazał dobre dopasowanie do danych historycznych, jednak prognozy należy traktować z ostrożnością z powodu ograniczonej liczby danych.

Lessons Learned / doubts:

  1. Większa ilość obserwacji mogłaby być bardziej efektywna przy analizie i forecast’ingu szeregu czasowego.

  2. Mała ilość anomalii oraz stabilne składowe (trend i sezonowość) powodują, że model przewiduje stabilną kontunaucję wartości, zgodnie ze składowymi

  3. Model ARIMA może teoretycznie przyjąć zbiór, który jest “nie do końca stacjonarny”. Przed różnicowaniem, tylko jeden z testów (ADF) wskazał stacjonarność szeregu. Dlatego wątpliwym jest, który z szeregów (zróżniocowany lub nie) powinien zostać użyty przy predykcji przyszłych wartości.

Poniżej, dla referencji - predykcja wartości z modelem auto.arima, na podstawie szeregu NIE-zróżnicowanego.

model_not_diff <- auto.arima(ts)
summary(model_not_diff)
## Series: ts 
## ARIMA(2,1,1)(0,1,0)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.5960  0.2143  -0.9819
## s.e.  0.0888  0.0880   0.0292
## 
## sigma^2 = 132.3:  log likelihood = -504.92
## AIC=1017.85   AICc=1018.17   BIC=1029.35
## 
## Training set error measures:
##                  ME     RMSE     MAE      MPE     MAPE     MASE        ACF1
## Training set 1.3423 10.84619 7.86754 0.420698 2.800458 0.245628 -0.00124847
forecasted_not_diff <- forecast(model_not_diff, h = 12)
plot(forecasted_not_diff)

Wykres przewidywanych wartości pozwala wysunąć bardzo podobne wnioski do predykcji na szeregu zróżnicowanym (stacjonarnym). Jak widać, model auto.arima sam dobrał parametr i(d) = 1, więc sam zróżnicował szereg.