Prognozowanie i symulacje - raport z projektu
Bukowska, Fleks, Mogielnicki, Pek
2024-06-15
Informacje o raporcie
Cel raportu
Celem naszego raportu jest wykonanie prognozy wielkości zarejestrowanych bezrobotnych na Pomorzu w 2023 roku na podstawie wybranych zmiennych. Zakresem danych każdej zmiennej w naszym zbiorze danych jest 2004 do 2022 roku. Dane do naszego zbioru danych zostały pobrane ze strony https://bdl.stat.gov.pl/bdl/dane/podgrup/temat.
Wyjaśnienie zmiennych
Bezrobotni_zarejestrowani (zmienna objaśniająca) - bezrobotni zarejestrowani na Pomorzu
Okres - zakres lat między pierwszą, a ostatnią obserwacją w zbiorze (2004-2022),
Absolwenci_uczelni [10 tys.] - absolwenci uczelni na Pomorzu, wybraliśmy tą zmienną ponieważ powinno to wpływać na bardziej wykwalifikowaną siłę roboczą a więc również na poziom bezrobocia
Oferty_pracy - oferty pracy na Pomorzu, zmienna została wybrana ponieważ ilość ofert pracy jest zagadnieniem bliskim dla bezrobocia
Prz_mie_wydatki_na_1_osobę - przeciętne miesięczne wydatki na 1 osobę na Pomorzu, wybór tej zmiennej uzasadniliśmy powiązaniem bezrobocia z poziomem życia, a więc z m. in. przeciętnymi wydatkami
Prz_mie_dochód_rozp_na_1_osobę - przeciętny miesięczny dochód rozporządzalny na 1 osobę na Pomorzu, wybrany z tego samego powodu co przeciętne wydatki
Prz_mie_wyn_brutto - przeciętne miesięczne wynagrodzenie brutto na Pomorzu, wybrane z tego samego powodu co dwie poprzednie zmienne
Stwierdzone_choroby_zawodowe - stwierdzone choroby zawodowe na Pomorzu, zmienna wybrana ponieważ choroby zawodowe mogą wpływać na zdolności pracowników do pracy
Inflacja - odnotowany poziom inflacji na Pomorzu
Wydatki_na_Oświata_i_wychowanie - wydatki na oświatę i wychowanie na Pomorzu, w celu zbadania czy takie wydatki wpływają pozytywnie na stopę bezrobocia
Analiza istotności zmiennych objaśniających
Aby wykonać wspomnianą prognozę wpierw utworzyliśmy pierwszy podstawowy model, aby określić nasze dalsze kroki. Przy tworzeniu modelu oraz weryfikacji hipotez statystycznych zostanie zastotsowany poziom istotności wynoszący 0,05.
Wstępny model
# wstepny model
model1 <- lm(Bezrobotni_zarejestrowani ~ Absolwenci_uczelni + Oferty_pracy + Prz_mie_wydatki_na_1_osobę + Prz_mie_dochód_rozp_na_1_osobę + Prz_mie_wyn_brutto + Stwierdzone_choroby_zawodowe + Inflacja + Wydatki_na_Oświata_i_wychowanie, data=df)
summary(model1)##
## Call:
## lm(formula = Bezrobotni_zarejestrowani ~ Absolwenci_uczelni +
## Oferty_pracy + Prz_mie_wydatki_na_1_osobę + Prz_mie_dochód_rozp_na_1_osobę +
## Prz_mie_wyn_brutto + Stwierdzone_choroby_zawodowe + Inflacja +
## Wydatki_na_Oświata_i_wychowanie, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37520 -6848 -1068 7009 21025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.253e+05 2.113e+05 2.013 0.07181 .
## Absolwenci_uczelni 3.356e+02 4.275e+02 0.785 0.45055
## Oferty_pracy -1.601e+01 4.629e+00 -3.458 0.00614 **
## Prz_mie_wydatki_na_1_osobę -7.040e+01 1.726e+02 -0.408 0.69197
## Prz_mie_dochód_rozp_na_1_osobę -1.122e+02 9.218e+01 -1.217 0.25147
## Prz_mie_wyn_brutto 4.470e+01 2.710e+01 1.649 0.13015
## Stwierdzone_choroby_zawodowe 1.428e+02 2.831e+02 0.504 0.62493
## Inflacja -2.222e+03 2.180e+03 -1.019 0.33211
## Wydatki_na_Oświata_i_wychowanie -2.130e-03 2.232e-03 -0.954 0.36241
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16840 on 10 degrees of freedom
## Multiple R-squared: 0.9003, Adjusted R-squared: 0.8205
## F-statistic: 11.28 on 8 and 10 DF, p-value: 0.0004264Zlogarytmowanie wybranych zmiennych
Uznaliśmy następnie, że wszystkie nasze zmienne należy zlogarytmować, z uwagi na to, że pierwotna postać modelu liniowego wskazuje na istotność statystyczną tylko jednej ze zmiennych, a model opracowany przy wykorzystaniu logarytmów często cechuje się większą dokładnością oraz może pozwolić na uwzględnienie w modelu większej ilości zmiennych.
variable_to_log <- c("Absolwenci_uczelni", "Bezrobotni_zarejestrowani", "Oferty_pracy", "Prz_mie_wydatki_na_1_osobę", "Prz_mie_dochód_rozp_na_1_osobę", "Prz_mie_wyn_brutto", "Stwierdzone_choroby_zawodowe", "Inflacja", "Wydatki_na_Oświata_i_wychowanie")
df_log <- df
df_log[variable_to_log] <- lapply(df[variable_to_log], log)Model ze zlogarytmowanymi zmiennymi
# model ze zlogarytmowanymi zmiennymi
model_log <- lm(Bezrobotni_zarejestrowani ~ Absolwenci_uczelni + Oferty_pracy + Prz_mie_wydatki_na_1_osobę + Prz_mie_dochód_rozp_na_1_osobę + Prz_mie_wyn_brutto + Stwierdzone_choroby_zawodowe + Inflacja + Wydatki_na_Oświata_i_wychowanie, data=df_log)
summary(model_log)##
## Call:
## lm(formula = Bezrobotni_zarejestrowani ~ Absolwenci_uczelni +
## Oferty_pracy + Prz_mie_wydatki_na_1_osobę + Prz_mie_dochód_rozp_na_1_osobę +
## Prz_mie_wyn_brutto + Stwierdzone_choroby_zawodowe + Inflacja +
## Wydatki_na_Oświata_i_wychowanie, data = df_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31870 -0.04605 -0.00378 0.07283 0.15017
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.98259 7.96076 3.138 0.01054 *
## Absolwenci_uczelni 0.86478 0.40450 2.138 0.05824 .
## Oferty_pracy -0.33427 0.09378 -3.565 0.00514 **
## Prz_mie_wydatki_na_1_osobę -0.61678 1.69031 -0.365 0.72279
## Prz_mie_dochód_rozp_na_1_osobę -1.83582 1.32180 -1.389 0.19502
## Prz_mie_wyn_brutto 0.96659 0.96952 0.997 0.34229
## Stwierdzone_choroby_zawodowe -0.11106 0.24347 -0.456 0.65803
## Inflacja -4.08901 1.94879 -2.098 0.06226 .
## Wydatki_na_Oświata_i_wychowanie 0.80480 0.49150 1.637 0.13258
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.146 on 10 degrees of freedom
## Multiple R-squared: 0.9419, Adjusted R-squared: 0.8955
## F-statistic: 20.28 on 8 and 10 DF, p-value: 3.186e-05Usuwanie nieistotnych zmiennych z modelu
W pierwotnej postaci modelu tylko jedna ze zmiennych objaśniających cechowała się istotnością statystyczną przy przyjętym poziomie istotności. Z tego powodu z naszego modelu utworzonego na podstawie zmiennych zlogarytmowanych usuwaliśmy po kolei zmienne nieistotne których wartośc p była największa, aby określić, które z nich mają istotny wpływ na zmienną objaśnianą.
model_log <- update(model_log, . ~ . -Prz_mie_wydatki_na_1_osobę)
summary(model_log)##
## Call:
## lm(formula = Bezrobotni_zarejestrowani ~ Absolwenci_uczelni +
## Oferty_pracy + Prz_mie_dochód_rozp_na_1_osobę + Prz_mie_wyn_brutto +
## Stwierdzone_choroby_zawodowe + Inflacja + Wydatki_na_Oświata_i_wychowanie,
## data = df_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31731 -0.04827 -0.00187 0.07223 0.14672
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.22777 7.61339 3.314 0.00691 **
## Absolwenci_uczelni 0.86589 0.38822 2.230 0.04749 *
## Oferty_pracy -0.32331 0.08526 -3.792 0.00299 **
## Prz_mie_dochód_rozp_na_1_osobę -2.20445 0.81810 -2.695 0.02086 *
## Prz_mie_wyn_brutto 0.85551 0.88347 0.968 0.35368
## Stwierdzone_choroby_zawodowe -0.09893 0.23149 -0.427 0.67736
## Inflacja -4.33669 1.75330 -2.473 0.03093 *
## Wydatki_na_Oświata_i_wychowanie 0.80643 0.47172 1.710 0.11537
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1401 on 11 degrees of freedom
## Multiple R-squared: 0.9412, Adjusted R-squared: 0.9037
## F-statistic: 25.14 on 7 and 11 DF, p-value: 6.345e-06model_log <- update(model_log, . ~ . -Stwierdzone_choroby_zawodowe)
summary(model_log)##
## Call:
## lm(formula = Bezrobotni_zarejestrowani ~ Absolwenci_uczelni +
## Oferty_pracy + Prz_mie_dochód_rozp_na_1_osobę + Prz_mie_wyn_brutto +
## Inflacja + Wydatki_na_Oświata_i_wychowanie, data = df_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.308483 -0.041041 -0.007914 0.070396 0.158781
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.41924 7.11899 3.430 0.00498 **
## Absolwenci_uczelni 0.82566 0.36358 2.271 0.04237 *
## Oferty_pracy -0.31974 0.08191 -3.903 0.00210 **
## Prz_mie_dochód_rozp_na_1_osobę -2.24632 0.78407 -2.865 0.01422 *
## Prz_mie_wyn_brutto 0.99008 0.79684 1.243 0.23778
## Inflacja -4.34530 1.69243 -2.567 0.02466 *
## Wydatki_na_Oświata_i_wychowanie 0.79142 0.45411 1.743 0.10691
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1352 on 12 degrees of freedom
## Multiple R-squared: 0.9402, Adjusted R-squared: 0.9103
## F-statistic: 31.44 on 6 and 12 DF, p-value: 1.154e-06model_log <- update(model_log, . ~ . -Prz_mie_wyn_brutto)
summary(model_log)##
## Call:
## lm(formula = Bezrobotni_zarejestrowani ~ Absolwenci_uczelni +
## Oferty_pracy + Prz_mie_dochód_rozp_na_1_osobę + Inflacja +
## Wydatki_na_Oświata_i_wychowanie, data = df_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31417 -0.04048 0.03532 0.07321 0.14908
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.73039 7.13270 3.187 0.007147 **
## Absolwenci_uczelni 0.52931 0.28010 1.890 0.081301 .
## Oferty_pracy -0.33118 0.08308 -3.986 0.001551 **
## Prz_mie_dochód_rozp_na_1_osobę -1.33334 0.27925 -4.775 0.000363 ***
## Inflacja -3.78382 1.66475 -2.273 0.040649 *
## Wydatki_na_Oświata_i_wychowanie 0.91611 0.45205 2.027 0.063727 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.138 on 13 degrees of freedom
## Multiple R-squared: 0.9325, Adjusted R-squared: 0.9065
## F-statistic: 35.91 on 5 and 13 DF, p-value: 3.642e-07model_log <- update(model_log, . ~ . -Absolwenci_uczelni)
summary(model_log)##
## Call:
## lm(formula = Bezrobotni_zarejestrowani ~ Oferty_pracy + Prz_mie_dochód_rozp_na_1_osobę +
## Inflacja + Wydatki_na_Oświata_i_wychowanie, data = df_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31137 -0.07061 0.04803 0.08922 0.16263
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.3169 6.7705 4.330 0.000692 ***
## Oferty_pracy -0.3811 0.0857 -4.447 0.000553 ***
## Prz_mie_dochód_rozp_na_1_osobę -1.3048 0.3034 -4.301 0.000732 ***
## Inflacja -5.2890 1.5904 -3.326 0.005001 **
## Wydatki_na_Oświata_i_wychowanie 1.0902 0.4815 2.264 0.039964 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1502 on 14 degrees of freedom
## Multiple R-squared: 0.9139, Adjusted R-squared: 0.8894
## F-statistic: 37.17 on 4 and 14 DF, p-value: 2.585e-07W finalnej wersji modelu objaśniającego kształtowanie się liczby zarejestrowanych osób bezrobotnych na Pomorzu znalazły się następujące zmienne objaśniające: -Oferty_pracy -Prz_mie_dochód_rozp_na_1_osobę -Inflacja -Wydatki_na_Oświata_i_wychowanie
Wyliczenie błędów ex-ante
# utworzenie funkcji wyliczajacej blad RMSPE
fun_rmspe <- function(actual, predicted) {
rmse_value <- sqrt(mean((actual - predicted)^2))
rmspe_value <- (rmse_value / mean(actual)) * 100
}
# obliczanie bledow
MAE <- mean(abs(df_log$Bezrobotni_zarejestrowani - model_log$fitted.values))
MSE <- mean((df_log$Bezrobotni_zarejestrowani - model_log$fitted.values)^2)
RMSE <- sqrt(MSE)
MAPE <- mean(abs((df_log$Bezrobotni_zarejestrowani - model_log$fitted.values) / df_log$Bezrobotni_zarejestrowani)) * 100
rmspe_value <- fun_rmspe(df_log$Bezrobotni_zarejestrowani, model_log$fitted.values)
# wyświetlenie bledow
cat("Średni absolutny błąd (MAE):", MAE, "\n")## Średni absolutny błąd (MAE): 0.1040787cat("Błąd średniokwadratowy (MSE):", MSE, "\n")## Błąd średniokwadratowy (MSE): 0.01661948cat("Pierwiastek z błędu średniokwadratowego (RMSE):", RMSE, "\n")## Pierwiastek z błędu średniokwadratowego (RMSE): 0.1289166cat("Średni absolutny procentowy błąd (MAPE):", MAPE, "%\n")## Średni absolutny procentowy błąd (MAPE): 0.934299 %cat("Pierwiastek z błędu średniokwadratowego procentowego (RMSPE): ", rmspe_value, "\n")## Pierwiastek z błędu średniokwadratowego procentowego (RMSPE): 1.141211Błąd RMSPE dla przewidywanych wartości zmiennej objaśnianej (Bezrobotni_zarejestrowani) wyniósł około 1,14. Oznacza to, że wartości przewidywane od rzeczywistych różnią się o 1,14%, co wskazuje na wysoką dokładność modelu. Model prognozuje zmienną objaśnianą z bardzo małym średnim procentowym błędem, co oznacza, że prognozowane wartości są bardzo bliskie rzeczywistym wartościom.
Testy statystyczne
Test na nieliniowość (kwadraty):
Hipoteza zerowa dla testu na nieliniowość kwadratów mówi, że w danych nie występuje nieliniowość, która może zostać wykryta przy dodaniu kwadratów do modelu. Natomiast hipoteza alternatywna mówi, że dodanie kwadratów do modelu istotnie wpłynie na detekcję nieliniowości dotyczącą modelu. Wartość p-value, wynosząca około 0,07, jest nieznacznie wyższa od przyjętego poziomu istotności wynoszącego około 0,05, co nie pozwala na odrzucenie hipotezy zerowej.
model_log_test_kwadraty <- lm(Bezrobotni_zarejestrowani ~ Oferty_pracy + Prz_mie_dochód_rozp_na_1_osobę + Inflacja + Wydatki_na_Oświata_i_wychowanie + I(Oferty_pracy^2) + I(Prz_mie_dochód_rozp_na_1_osobę^2) + I(Inflacja^2) + I(Wydatki_na_Oświata_i_wychowanie^2) , data = df_log)
anova(model_log, model_log_test_kwadraty)## Analysis of Variance Table
##
## Model 1: Bezrobotni_zarejestrowani ~ Oferty_pracy + Prz_mie_dochód_rozp_na_1_osobę +
## Inflacja + Wydatki_na_Oświata_i_wychowanie
## Model 2: Bezrobotni_zarejestrowani ~ Oferty_pracy + Prz_mie_dochód_rozp_na_1_osobę +
## Inflacja + Wydatki_na_Oświata_i_wychowanie + I(Oferty_pracy^2) +
## I(Prz_mie_dochód_rozp_na_1_osobę^2) + I(Inflacja^2) + I(Wydatki_na_Oświata_i_wychowanie^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 14 0.31577
## 2 10 0.14304 4 0.17273 3.0189 0.07124 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test na nieliniowość (logarytmy):
model_log_test_log <- lm(Bezrobotni_zarejestrowani ~ log(Oferty_pracy) + log(Prz_mie_dochód_rozp_na_1_osobę) + log(Inflacja) + log(Wydatki_na_Oświata_i_wychowanie) , data = df_log)
anova(model_log, model_log_test_log)## Analysis of Variance Table
##
## Model 1: Bezrobotni_zarejestrowani ~ Oferty_pracy + Prz_mie_dochód_rozp_na_1_osobę +
## Inflacja + Wydatki_na_Oświata_i_wychowanie
## Model 2: Bezrobotni_zarejestrowani ~ log(Oferty_pracy) + log(Prz_mie_dochód_rozp_na_1_osobę) +
## log(Inflacja) + log(Wydatki_na_Oświata_i_wychowanie)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 14 0.31577
## 2 14 0.31698 0 -0.001212Wybrane testy statystyczne
Poniżej zostały przedstawione wyniki różnych testów statystycznych stosowanych w celu badania różnych właściwości statystycznych modeli.
dwtest(model_log)##
## Durbin-Watson test
##
## data: model_log
## DW = 1.4148, p-value = 0.01454
## alternative hypothesis: true autocorrelation is greater than 0bptest(model_log)##
## studentized Breusch-Pagan test
##
## data: model_log
## BP = 1.2352, df = 4, p-value = 0.8723bgtest(model_log)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_log
## LM test = 1.6787, df = 1, p-value = 0.1951resettest(model_log)##
## RESET test
##
## data: model_log
## RESET = 1.5114, df1 = 2, df2 = 12, p-value = 0.2598shapiro.test(resid(model_log))##
## Shapiro-Wilk normality test
##
## data: resid(model_log)
## W = 0.90795, p-value = 0.06789jarque.bera.test(resid(model_log))##
## Jarque Bera Test
##
## data: resid(model_log)
## X-squared = 2.8455, df = 2, p-value = 0.2411durbinWatsonTest(model_log)## lag Autocorrelation D-W Statistic p-value
## 1 0.2722291 1.414784 0.042
## Alternative hypothesis: rho != 0Otrzymane wyniki należy interpretować następująco:
Test Durbina Watsona służy do weryfikacji występowania autokorelacji reszt. Jego hipoteza zerowa mówi, że autokorelacja reszt nie występuje, natomiast hipoteza alternatywna wskazuje na jej występowanie. Wartość p-value dla tego modelu, wynosząca około 0.014 oraz niższa od przyjętego poziomu istotności, nie pozwala na odrzucenie hipotezy alternatywnej.
Test Breuscha-Pagana służy do wykrywania tego, czy w modelu dochodzi do heteroskedastyczności, czyli sytuacji, w której wariancja reszt dotycząca modelu nie jest stała. Hipoteza zerowa tego testu mówi, że heteroskedastycznośc reszt nie występuje, a więc występuje homoskedastyczność. Na podstawie wartości p-value dla tego testu, wynoszącej 0.8723, nie można odrzucić hipotezy zerowej.
Test Breuscha-Godfreya jest kolejnym testem służącym do wykrywania występowania autokorelacji reszt, jednak w przeciwieństwie do testu Durbina Watsona badającego wyłącznie autokorelację pierwszego rzędu, ten test bada również występowanie autokorelacji wyższych rzędów. Jego hipoteza zerowa, której nie można odrzucić na podstawie otrzymanej wartości p-value nie można odrzucić, mówi o braku występowania autokorelacji w analizowanym zakresie.
Test RESET służy do sprawdzenia poprawności specyfikacji modelu. Jego hipoteza zerowa, której przy przyjętym poziomie istotności nie można odrzucić, mowi, że specyfikacja modelu jest poprawna.
Test Shapiro-Wilka sprawdza, czy wykorzystywane dane cechują się rozkładem zgodnym z rozkładem normalnym. Hipoteza zerowa mówi, że rozkład danych jest normalny. Przy przyjętym poziomie istotności nie można jej odrzucić, ale wartośc p-value wynosząca około 0.07 jest nieznacznie wyższa od przyjętego poziomu istotności wynoszącego 0.05, co można uzasadnić tym tym, że wykorzystywane dane są wartościami zlogarytmowanymi.
Test Jarque Bery, podobnie jak test Shapiro Wilka, weryfikuje normalność rozkładu dotyczącą wykorzystanych danych, ale jest często stosowany dla szerszego zakresu danych. Jego hipotezy zerowej, zgodnej z hipotezą zerową dla testu Shapiro Wilka oraz mówiącej, że wykorzystywana próbka jest zgodna z rozkłądem normalnym, nie można odrzucić.
Wizualizacja wykresów zmiennych
# utworzenie zbioru danych, który zawiera tylko okres i zmienne, które okazały się być w modelu istotne statystycznie
df_ost <- subset(df_log, select = c(Bezrobotni_zarejestrowani, Oferty_pracy, Prz_mie_dochód_rozp_na_1_osobę, Inflacja, Wydatki_na_Oświata_i_wychowanie))
lista_wykresow <- list() # pusta lista do przechowywania wykresów
# tworzenie wykresów dla zmiennej objaśnianej i zmiennych objaśniających, które zostały uznane za istotne statystycznie w naszym modelu
for (colname in colnames(df_ost)) {
nazwa_wykresu <- paste("wykres_", colname, sep = "")
# Tworzenie wykresu i zapisanie go do listy
assign(nazwa_wykresu,
ggplot(df_ost, aes_string(x = df_log$Okres, y = colname)) +
geom_line() +
geom_point() +
labs(title = paste(colname, "w czasie"), x = "Rok", y = colname) +
theme_minimal()
)
lista_wykresow[[nazwa_wykresu]] <- get(nazwa_wykresu) #odanie utworzonego wykresu do listy z wykresami
}## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# wyświetlenie wszystkich wykresów z listy
for (wykres in lista_wykresow) {
print(wykres)
}
Korelacja
Kolejnym krokiem będzie utworzenie macierzy korelacji dotyczącej wykorzystywanych zmiennych.
correlation <- data.frame(cor(df_ost))
correlation## Bezrobotni_zarejestrowani Oferty_pracy
## Bezrobotni_zarejestrowani 1.0000000 -0.78432858
## Oferty_pracy -0.7843286 1.00000000
## Prz_mie_dochód_rozp_na_1_osobę -0.8728012 0.65130527
## Inflacja -0.3340565 -0.08811464
## Wydatki_na_Oświata_i_wychowanie -0.7223021 0.51645926
## Prz_mie_dochód_rozp_na_1_osobę Inflacja
## Bezrobotni_zarejestrowani -0.8728012 -0.33405651
## Oferty_pracy 0.6513053 -0.08811464
## Prz_mie_dochód_rozp_na_1_osobę 1.0000000 0.30735531
## Inflacja 0.3073553 1.00000000
## Wydatki_na_Oświata_i_wychowanie 0.8830197 0.49873667
## Wydatki_na_Oświata_i_wychowanie
## Bezrobotni_zarejestrowani -0.7223021
## Oferty_pracy 0.5164593
## Prz_mie_dochód_rozp_na_1_osobę 0.8830197
## Inflacja 0.4987367
## Wydatki_na_Oświata_i_wychowanie 1.0000000Współliniowość
vif <- vif(lm(Bezrobotni_zarejestrowani ~ ., data = df_ost))
vif## Oferty_pracy Prz_mie_dochód_rozp_na_1_osobę
## 2.070985 6.045038
## Inflacja Wydatki_na_Oświata_i_wychowanie
## 1.729185 6.141973Prognozowanie zmiennych objaśniających
Utworzenie z istotnych statystycznie zmiennych szeregi czasowe
Indeks = 1:19
df_log <- df_log %>%
mutate(Indeks) %>%
select(Indeks, everything())
ym("2004-01")## [1] "2004-01-01"# tworzenie zmiennych
log_oferty_pracy <- df_log$Oferty_pracy
log_prz_mie_dochód_rozp_na_1_osobę <- df_log$Prz_mie_dochód_rozp_na_1_osobę
log_inflacja <- df_log$Inflacja
log_wydatki_na_Oświata_i_wychowanie <- df_log$Wydatki_na_Oświata_i_wychowanie
# zamienienie tych zmiennych na szeregi czasowe
log_oferty_pracy.ts <- ts(log_oferty_pracy, frequency=1, start=c(2004))
log_prz_mie_dochód_rozp_na_1_osobę.ts <- ts(log_prz_mie_dochód_rozp_na_1_osobę, frequency=1, start=c(2004))
log_inflacja.ts <- ts(log_inflacja, frequency=1, start=c(2004))
log_wydatki_na_Oświata_i_wychowanie.ts <- ts(log_wydatki_na_Oświata_i_wychowanie, frequency=1, start=c(2004))
# utworzenie zbioru danych z utworzonych szeregow czasowych
df_ts <- cbind(log_oferty_pracy.ts, log_prz_mie_dochód_rozp_na_1_osobę.ts, log_inflacja.ts, log_wydatki_na_Oświata_i_wychowanie.ts)
okres_prognozy <- data.frame(Indeks = seq(20,20))W celu wykonania prognozy dla liczby zarejestrowanych osób bezrobotnych na Pomorzu w roku 2023 najpierw należało wykonać prognozy dla zmiennych objaśniających które znalazły się w naszym finalnym modelu.
Prognozowanie - trend liniowy
rmspe_values_lm <- list() #lista na wyniki RMSPE
lista_wykresow_lm <- list()
for (var in colnames(df_ts)) {
# utworzenie modeli
print("Model")
print(var)
model_progn <- lm(df_ts[, var] ~ Indeks, data = df_ts)
print(summary(model_progn))
# przeprowadzenie testow statystycznych
print("Testy statystyczne")
print(dwtest(model_progn))
print(bgtest(model_progn))
print(bptest(model_progn))
print(resettest(model_progn))
# zaprognozowanie zmiennych i okreslenie dokladnosci prognozy
print("Dokładność prognozy")
var_prognoza <- forecast(model_progn, h=1, newdata=okres_prognozy)
errors <- accuracy(var_prognoza) #bledy prognozy
print(errors)
rmspe_value <- fun_rmspe(df_ts[, var], var_prognoza$fitted)
cat("RMSPE: ", rmspe_value, "\n")
rmspe_values_lm[[var]] <- paste("RMSPE (lm): ", rmspe_value)
cat("-------------------------------------------------------------------------------------------------\n")
# utworzenie dataframe z danymi potrzebnymi do stworzenia wykresow
wykres_dane_lm <- data.frame(okres = c(df_log$Okres, 2023), actual = c(df_ts[, var], NA), predicted = c(var_prognoza$fitted, var_prognoza$mean[1]))
# narysowanie wykresow
nazwa_wykresu <- paste("wykres_lm_", var, sep = "")
assign(nazwa_wykresu,
ggplot(wykres_dane_lm, aes(x = okres)) +
geom_line(aes(y = actual), color = "blue", linetype = "solid", size = 1) +
geom_line(aes(y = predicted), color = "red", linetype = "dashed", size = 1) +
labs(title = paste("Rzeczywiste vs przewidziane wartości trendem liniowym dla:", var),
x = "Rok",
y = var) +
theme_minimal()
)
lista_wykresow_lm[[nazwa_wykresu]] <- get(nazwa_wykresu)
}## [1] "Model"
## [1] "log_oferty_pracy.ts"
##
## Call:
## lm(formula = df_ts[, var] ~ Indeks, data = df_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.93243 -0.26843 -0.04311 0.26539 1.05134
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.98307 0.22479 31.064 < 2e-16 ***
## Indeks 0.06746 0.01972 3.422 0.00325 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4707 on 17 degrees of freedom
## Multiple R-squared: 0.4078, Adjusted R-squared: 0.373
## F-statistic: 11.71 on 1 and 17 DF, p-value: 0.003252
##
## [1] "Testy statystyczne"
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.88294, p-value = 0.001204
## alternative hypothesis: true autocorrelation is greater than 0
##
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 3.5624, df = 1, p-value = 0.0591
##
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 2.82, df = 1, p-value = 0.09309
##
##
## RESET test
##
## data: model_progn
## RESET = 0.59061, df1 = 2, df2 = 15, p-value = 0.5664
##
## [1] "Dokładność prognozy"
## ME RMSE MAE MPE MAPE MASE
## Training set -1.402387e-16 0.4452449 0.3387036 -0.3637291 4.502558 0.7181763
## RMSPE: 5.814369
## -------------------------------------------------------------------------------------------------## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## [1] "Model"
## [1] "log_prz_mie_dochód_rozp_na_1_osobę.ts"
##
## Call:
## lm(formula = df_ts[, var] ~ Indeks, data = df_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12374 -0.04109 0.01349 0.06017 0.07889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.698524 0.031302 214.00 < 2e-16 ***
## Indeks 0.049710 0.002745 18.11 1.51e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06554 on 17 degrees of freedom
## Multiple R-squared: 0.9507, Adjusted R-squared: 0.9478
## F-statistic: 327.9 on 1 and 17 DF, p-value: 1.51e-12
##
## [1] "Testy statystyczne"
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.52691, p-value = 7.88e-06
## alternative hypothesis: true autocorrelation is greater than 0
##
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 7.5447, df = 1, p-value = 0.006019
##
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 5.6636, df = 1, p-value = 0.01732
##
##
## RESET test
##
## data: model_progn
## RESET = 17.123, df1 = 2, df2 = 15, p-value = 0.0001343
##
## [1] "Dokładność prognozy"
## ME RMSE MAE MPE MAPE MASE
## Training set 0 0.06199876 0.05375325 -0.008624115 0.7584353 0.2372324
## RMSPE: 0.8616177
## -------------------------------------------------------------------------------------------------
## [1] "Model"
## [1] "log_inflacja.ts"
##
## Call:
## lm(formula = df_ts[, var] ~ Indeks, data = df_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.038078 -0.017008 0.003857 0.008033 0.088172
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.614761 0.013588 339.61 <2e-16 ***
## Indeks 0.001705 0.001192 1.43 0.171
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02845 on 17 degrees of freedom
## Multiple R-squared: 0.1074, Adjusted R-squared: 0.05491
## F-statistic: 2.046 on 1 and 17 DF, p-value: 0.1708
##
## [1] "Testy statystyczne"
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.69228, p-value = 0.000124
## alternative hypothesis: true autocorrelation is greater than 0
##
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 7.888, df = 1, p-value = 0.004977
##
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 3.35, df = 1, p-value = 0.0672
##
##
## RESET test
##
## data: model_progn
## RESET = 19.549, df1 = 2, df2 = 15, p-value = 6.635e-05
##
## [1] "Dokładność prognozy"
## ME RMSE MAE MPE MAPE MASE
## Training set 4.674659e-17 0.02691419 0.01892314 -0.003332819 0.4071567 1.067909
## RMSPE: 0.5810733
## -------------------------------------------------------------------------------------------------
## [1] "Model"
## [1] "log_wydatki_na_Oświata_i_wychowanie.ts"
##
## Call:
## lm(formula = df_ts[, var] ~ Indeks, data = df_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21971 -0.06207 0.04387 0.06931 0.12904
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.95459 0.04697 360.942 < 2e-16 ***
## Indeks 0.02756 0.00412 6.691 3.8e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09836 on 17 degrees of freedom
## Multiple R-squared: 0.7248, Adjusted R-squared: 0.7086
## F-statistic: 44.77 on 1 and 17 DF, p-value: 3.796e-06
##
## [1] "Testy statystyczne"
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.53881, p-value = 9.952e-06
## alternative hypothesis: true autocorrelation is greater than 0
##
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 6.1692, df = 1, p-value = 0.013
##
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 0.61807, df = 1, p-value = 0.4318
##
##
## RESET test
##
## data: model_progn
## RESET = 29.125, df1 = 2, df2 = 15, p-value = 6.833e-06
##
## [1] "Dokładność prognozy"
## ME RMSE MAE MPE MAPE MASE
## Training set 0 0.09303873 0.0778458 -0.002943661 0.4528919 0.6137119
## RMSPE: 0.5399735
## -------------------------------------------------------------------------------------------------rmspe_values_lm## $log_oferty_pracy.ts
## [1] "RMSPE (lm): 5.8143692677273"
##
## $log_prz_mie_dochód_rozp_na_1_osobę.ts
## [1] "RMSPE (lm): 0.861617734292072"
##
## $log_inflacja.ts
## [1] "RMSPE (lm): 0.581073268277916"
##
## $log_wydatki_na_Oświata_i_wychowanie.ts
## [1] "RMSPE (lm): 0.539973489648822"for (wykres in lista_wykresow_lm) {
print(wykres)
}## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
Prognozowanie - trend kwadratowy
rmspe_values_lm2 <- list() #lista na wyniki RMSPE
lista_wykresow_lm2 <- list()
for (var2 in colnames(df_ts)) {
# utworzenie modeli
print("Model")
print(var2)
model_progn <- lm(df_ts[, var2] ~ Indeks + I(Indeks^2), data = df_ts)
print(summary(model_progn))
# przeprowadzenie testow statystycznych
print("Testy statystyczne")
print(dwtest(model_progn))
print(bgtest(model_progn))
print(bptest(model_progn))
print(resettest(model_progn))
# zaprognozowanie zmiennych i okreslenie dokladnosci prognozy
print("Dokładność prognozy")
var2_prognoza <- forecast(model_progn, h=1, newdata=okres_prognozy)
errors <- accuracy(var2_prognoza) #bledy prognozy
print(errors)
rmspe_value <- fun_rmspe(df_ts[, var2], var2_prognoza$fitted)
cat("RMSPE: ", rmspe_value, "\n")
rmspe_values_lm2[[var2]] <- paste("RMSPE (lm): ", rmspe_value)
cat("-------------------------------------------------------------------------------------------------\n")
# utworzenie dataframe z danymi potrzebnymi do stworzenia wykresow
wykres_dane_lm2 <- data.frame(okres = c(df_log$Okres, 2023), actual = c(df_ts[, var2], NA), predicted = c(var2_prognoza$fitted, var2_prognoza$mean[1]))
# narysowanie wykresow
nazwa_wykresu <- paste("wykres_lm2_", var2, sep = "")
assign(nazwa_wykresu,
ggplot(wykres_dane_lm2, aes(x = okres)) +
geom_line(aes(y = actual), color = "blue", linetype = "solid", size = 1) +
geom_line(aes(y = predicted), color = "red", linetype = "dashed", size = 1) +
labs(title = paste("Rzeczywiste vs przewidziane wartości trendem kwadratowym wartości dla:", var2),
x = "Rok",
y = var2) +
theme_minimal()
)
lista_wykresow_lm2[[nazwa_wykresu]] <- get(nazwa_wykresu)
}## [1] "Model"
## [1] "log_oferty_pracy.ts"
##
## Call:
## lm(formula = df_ts[, var2] ~ Indeks + I(Indeks^2), data = df_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7040 -0.1919 -0.1058 0.2967 1.0782
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.669542 0.358746 18.591 2.94e-12 ***
## Indeks 0.157039 0.082601 1.901 0.0754 .
## I(Indeks^2) -0.004479 0.004012 -1.116 0.2808
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4673 on 16 degrees of freedom
## Multiple R-squared: 0.4506, Adjusted R-squared: 0.3819
## F-statistic: 6.561 on 2 and 16 DF, p-value: 0.0083
##
## [1] "Testy statystyczne"
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.89761, p-value = 0.0003738
## alternative hypothesis: true autocorrelation is greater than 0
##
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 4.3632, df = 1, p-value = 0.03672
##
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 3.0437, df = 2, p-value = 0.2183
##
##
## RESET test
##
## data: model_progn
## RESET = 10.331, df1 = 2, df2 = 14, p-value = 0.001754
##
## [1] "Dokładność prognozy"
## ME RMSE MAE MPE MAPE MASE
## Training set 9.349033e-17 0.4288583 0.3433828 -0.3327313 4.564371 0.728098
## RMSPE: 5.60038
## -------------------------------------------------------------------------------------------------
## [1] "Model"
## [1] "log_prz_mie_dochód_rozp_na_1_osobę.ts"
##
## Call:
## lm(formula = df_ts[, var2] ~ Indeks + I(Indeks^2), data = df_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.107123 -0.024431 0.007863 0.042299 0.051546
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.5838716 0.0367344 179.229 < 2e-16 ***
## Indeks 0.0824677 0.0084581 9.750 3.9e-08 ***
## I(Indeks^2) -0.0016379 0.0004109 -3.987 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04785 on 16 degrees of freedom
## Multiple R-squared: 0.9753, Adjusted R-squared: 0.9722
## F-statistic: 315.5 on 2 and 16 DF, p-value: 1.399e-13
##
## [1] "Testy statystyczne"
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.92763, p-value = 0.0005275
## alternative hypothesis: true autocorrelation is greater than 0
##
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 4.9038, df = 1, p-value = 0.0268
##
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 0.37069, df = 2, p-value = 0.8308
##
##
## RESET test
##
## data: model_progn
## RESET = 5.5751, df1 = 2, df2 = 14, p-value = 0.01656
##
## [1] "Dokładność prognozy"
## ME RMSE MAE MPE MAPE MASE
## Training set 0 0.04391371 0.03658238 -0.003939899 0.511205 0.1614512
## RMSPE: 0.6102837
## -------------------------------------------------------------------------------------------------
## [1] "Model"
## [1] "log_inflacja.ts"
##
## Call:
## lm(formula = df_ts[, var2] ~ Indeks + I(Indeks^2), data = df_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.02269 -0.01633 -0.00999 0.01438 0.05798
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.6562024 0.0182147 255.629 <2e-16 ***
## Indeks -0.0101357 0.0041939 -2.417 0.0280 *
## I(Indeks^2) 0.0005920 0.0002037 2.906 0.0103 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02373 on 16 degrees of freedom
## Multiple R-squared: 0.4158, Adjusted R-squared: 0.3427
## F-statistic: 5.693 on 2 and 16 DF, p-value: 0.01357
##
## [1] "Testy statystyczne"
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.86242, p-value = 0.0002442
## alternative hypothesis: true autocorrelation is greater than 0
##
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 6.6243, df = 1, p-value = 0.01006
##
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 4.4924, df = 2, p-value = 0.1058
##
##
## RESET test
##
## data: model_progn
## RESET = 16.903, df1 = 2, df2 = 14, p-value = 0.0001847
##
## [1] "Dokładność prognozy"
## ME RMSE MAE MPE MAPE MASE
## Training set 9.349282e-17 0.02177449 0.01797118 -0.002181217 0.3868959 1.014186
## RMSPE: 0.4701078
## -------------------------------------------------------------------------------------------------
## [1] "Model"
## [1] "log_wydatki_na_Oświata_i_wychowanie.ts"
##
## Call:
## lm(formula = df_ts[, var2] ~ Indeks + I(Indeks^2), data = df_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.20598 -0.06307 0.04180 0.06587 0.14277
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.9357441 0.0775951 218.258 <2e-16 ***
## Indeks 0.0329498 0.0178662 1.844 0.0837 .
## I(Indeks^2) -0.0002692 0.0008679 -0.310 0.7604
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1011 on 16 degrees of freedom
## Multiple R-squared: 0.7264, Adjusted R-squared: 0.6922
## F-statistic: 21.24 on 2 and 16 DF, p-value: 3.139e-05
##
## [1] "Testy statystyczne"
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.53723, p-value = 7.667e-07
## alternative hypothesis: true autocorrelation is greater than 0
##
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 6.8912, df = 1, p-value = 0.008662
##
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 2.7007, df = 2, p-value = 0.2592
##
##
## RESET test
##
## data: model_progn
## RESET = 27.033, df1 = 2, df2 = 14, p-value = 1.557e-05
##
## [1] "Dokładność prognozy"
## ME RMSE MAE MPE MAPE MASE
## Training set 0 0.09276014 0.07742067 -0.002919753 0.4501926 0.6103603
## RMSPE: 0.5383566
## -------------------------------------------------------------------------------------------------rmspe_values_lm2## $log_oferty_pracy.ts
## [1] "RMSPE (lm): 5.60037960941203"
##
## $log_prz_mie_dochód_rozp_na_1_osobę.ts
## [1] "RMSPE (lm): 0.610283665221309"
##
## $log_inflacja.ts
## [1] "RMSPE (lm): 0.470107826683551"
##
## $log_wydatki_na_Oświata_i_wychowanie.ts
## [1] "RMSPE (lm): 0.538356643913684"for (wykres in lista_wykresow_lm2) {
print(wykres)
}## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
Prognozowanie - Ruchoma zmienna ważona
rmspe_values_ma <- list() #lista na wyniki RMSPE
lista_wykresow_ma <- list()
#wykres_dane_ma <- data.frame(okres = df_log$Okres, actual = df_ost[[var]], predicted = ma_prognoza$fitted)
for (var in colnames(df_ost[-1])) {
# prognoza
#ma_model <- ma(df_ost[, var], order = 5)
ma_prognoza <- forecast(ma(df_ost[[var]], order = 3))
# wypisanie wynikow
print(var)
print(summary(ma_prognoza))
rmspe_value <- fun_rmspe(df_ost[[var]][3:19], ma_prognoza$fitted)
cat("RMSPE: ", rmspe_value, "\n")
rmspe_values_ma[[var]] <- paste("RMSPE (MA): ", rmspe_value)
cat("-------------------------------------------------------------------------------------------------\n")
# utworzenie dataframe z danymi potrzebnymi do stworzenia wykresów
ma_prognoza$fitted <- c(NA, NA, ma_prognoza$fitted)
wykres_dane_ma <- data.frame(okres = c(df_log$Okres, 2023), actual = c(df_ost[[var]], NA), predicted = c(ma_prognoza$fitted, ma_prognoza$mean[1]))
# narysowanie wykresu
nazwa_wykresu <- paste("wykres_ma_", var, sep = "")
assign(nazwa_wykresu,
ggplot(wykres_dane_ma, aes(x = okres)) +
geom_line(aes(y = actual), color = "blue", linetype = "solid", size = 1) +
geom_line(aes(y = predicted), color = "red", linetype = "dashed", size = 1) +
labs(title = paste("Rzeczywiste vs przewidziane wartości metodą MA dla:", var),
x = "Rok",
y = var) +
theme_minimal()
)
lista_wykresow_ma[[nazwa_wykresu]] <- get(nazwa_wykresu)
}## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series## [1] "Oferty_pracy"
##
## Forecast method: ETS(A,N,N)
##
## Model Information:
## ETS(A,N,N)
##
## Call:
## ets(y = object, lambda = lambda, biasadj = biasadj, allow.multiplicative.trend = allow.multiplicative.trend)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 6.8505
##
## sigma: 0.2777
##
## AIC AICc BIC
## 8.47553 10.32168 10.97517
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.06587926 0.2608533 0.1859605 0.8307823 2.407021 0.9412841
## ACF1
## Training set 0.4406536
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 19 7.970348 7.614462 8.326235 7.426067 8.514630
## 20 7.970348 7.467074 8.473622 7.200657 8.740040
## 21 7.970348 7.353976 8.586720 7.027689 8.913008
## 22 7.970348 7.258629 8.682068 6.881868 9.058829
## 23 7.970348 7.174626 8.766071 6.753396 9.187301
## 24 7.970348 7.098681 8.842016 6.637248 9.303448
## 25 7.970348 7.028842 8.911854 6.530439 9.410258
## 26 7.970348 6.963838 8.976859 6.431023 9.509673
## 27 7.970348 6.902784 9.037912 6.337650 9.603047
## 28 7.970348 6.845038 9.095658 6.249335 9.691362
## RMSPE: 6.905638
## -------------------------------------------------------------------------------------------------## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series## [1] "Prz_mie_dochód_rozp_na_1_osobę"
##
## Forecast method: ETS(M,A,N)
##
## Model Information:
## ETS(M,A,N)
##
## Call:
## ets(y = object, lambda = lambda, biasadj = biasadj, allow.multiplicative.trend = allow.multiplicative.trend)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.9975
##
## Initial states:
## l = 6.604
## b = 0.1039
##
## sigma: 0.0027
##
## AIC AICc BIC
## -80.17392 -74.71938 -76.00786
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.004500158 0.0173071 0.0140267 -0.06273178 0.1925447 0.2696514
## ACF1
## Training set 0.4259097
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 19 7.567982 7.541665 7.594299 7.527734 7.608230
## 20 7.595527 7.536698 7.654355 7.505556 7.685497
## 21 7.623071 7.524591 7.721552 7.472458 7.773684
## 22 7.650616 7.506362 7.794869 7.429999 7.871232
## 23 7.678160 7.482694 7.873626 7.379221 7.977100
## 24 7.705705 7.454077 7.957333 7.320873 8.090537
## 25 7.733250 7.420882 8.045617 7.255525 8.210975
## 26 7.760794 7.383404 8.138185 7.183625 8.337963
## 27 7.788339 7.341882 8.234796 7.105542 8.471136
## 28 7.815883 7.296516 8.335251 7.021579 8.610188
## RMSPE: 0.9574071
## -------------------------------------------------------------------------------------------------## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series## [1] "Inflacja"
##
## Forecast method: ETS(M,N,N)
##
## Model Information:
## ETS(M,N,N)
##
## Call:
## ets(y = object, lambda = lambda, biasadj = biasadj, allow.multiplicative.trend = allow.multiplicative.trend)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 4.6236
##
## sigma: 0.0026
##
## AIC AICc BIC
## -98.06280 -96.21665 -95.56316
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003062023 0.01138432 0.007771421 0.06555596 0.167614 0.9413895
## ACF1
## Training set 0.3948147
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 19 4.675641 4.659964 4.691318 4.651665 4.699616
## 20 4.675641 4.653471 4.697810 4.641736 4.709546
## 21 4.675641 4.648489 4.702792 4.634116 4.717165
## 22 4.675641 4.644289 4.706992 4.627693 4.723589
## 23 4.675641 4.640589 4.710693 4.622034 4.729248
## 24 4.675641 4.637243 4.714038 4.616917 4.734364
## 25 4.675641 4.634167 4.717114 4.612212 4.739069
## 26 4.675641 4.631303 4.719978 4.607833 4.743449
## 27 4.675641 4.628614 4.722668 4.603719 4.747562
## 28 4.675641 4.626070 4.725211 4.599829 4.751453
## RMSPE: 0.6136527
## -------------------------------------------------------------------------------------------------## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series## [1] "Wydatki_na_Oświata_i_wychowanie"
##
## Forecast method: ETS(A,N,N)
##
## Model Information:
## ETS(A,N,N)
##
## Call:
## ets(y = object, lambda = lambda, biasadj = biasadj, allow.multiplicative.trend = allow.multiplicative.trend)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 16.9603
##
## sigma: 0.0564
##
## AIC AICc BIC
## -45.71250 -43.86634 -43.21286
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03217074 0.05299475 0.04353641 0.1862162 0.2522055 0.9412576
## ACF1
## Training set 0.4781439
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 19 17.50711 17.43481 17.57941 17.39653 17.61768
## 20 17.50711 17.40486 17.60935 17.35074 17.66348
## 21 17.50711 17.38189 17.63233 17.31560 17.69862
## 22 17.50711 17.36252 17.65170 17.28597 17.72824
## 23 17.50711 17.34545 17.66877 17.25987 17.75434
## 24 17.50711 17.33002 17.68420 17.23628 17.77794
## 25 17.50711 17.31583 17.69838 17.21458 17.79964
## 26 17.50711 17.30263 17.71159 17.19438 17.81984
## 27 17.50711 17.29022 17.72399 17.17541 17.83881
## 28 17.50711 17.27849 17.73573 17.15747 17.85675
## RMSPE: 0.6259563
## -------------------------------------------------------------------------------------------------rmspe_values_ma## $Oferty_pracy
## [1] "RMSPE (MA): 6.90563841214671"
##
## $Prz_mie_dochód_rozp_na_1_osobę
## [1] "RMSPE (MA): 0.95740707253437"
##
## $Inflacja
## [1] "RMSPE (MA): 0.613652739293226"
##
## $Wydatki_na_Oświata_i_wychowanie
## [1] "RMSPE (MA): 0.625956299173227"for (wykres in lista_wykresow_ma) {
print(wykres)
}## Warning: Removed 1 row containing missing values (`geom_line()`).## Warning: Removed 2 rows containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Removed 2 rows containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Removed 2 rows containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Removed 2 rows containing missing values (`geom_line()`).
Prognozowanie - metoda Holta (pierwszy sposób)
rmspe_values_holt1 <- list() #lista na wyniki RMSPE
lista_wykresow_holt1 <- list()
for (var in colnames(df_ost[-1])) {
# prognoza
holt_model1 <- holt(df_ost[[var]], alpha=NULL, beta=NULL)
holt_prognoza1 <- forecast(holt_model1, h = 1)
# wypisanie wynikow
print("Model")
print(var)
print(summary(holt_prognoza1))
rmspe_value <- fun_rmspe(df_ost[[var]], holt_prognoza1$fitted)
cat("RMSPE: ", rmspe_value, "\n")
rmspe_values_holt1[[var]] <- paste("RMSPE (holt1): ", rmspe_value)
cat("-------------------------------------------------------------------------------\n")
# wyliczenie reszt
print("Testy dla reszt")
res <- holt_prognoza1$resid
plot(res, ylab = paste("Reszty (", var, ")", sep = "")) #wizualizacja reszt
# Testy dla reszt
print(acf(res, type = c("correlation"), plot = TRUE, ylab = paste("ACF (", var, ")", sep = "")))
print(pacf(res, plot = TRUE, ylab = paste("PACF (", var, ")", sep = "")))
print(adf.test(res))
print(Box.test (res, lag = , type = "Ljung"))
cat("-------------------------------------------------------------------------------\n")
# utworzenie dataframe z danymi potrzebnymi do stworzenia wykresow
wykres_dane_holt1 <- data.frame(okres = c(df_log$Okres, 2023), actual = c(df_ost[[var]], NA), predicted = c(holt_prognoza1$fitted, holt_prognoza1$mean))
wykres_prognozowane_dane_holt1 <- data.frame(Time = 2023,
Forecast = as.numeric(holt_prognoza1$mean),
Lo80 = holt_prognoza1$lower[, "80%"],
Hi80 = holt_prognoza1$upper[, "80%"],
Lo95 = holt_prognoza1$lower[, "95%"],
Hi95 = holt_prognoza1$upper[, "95%"])
# narysowanie wykresow
nazwa_wykresu <- paste("wykres_holt1_", var, sep = "")
assign(nazwa_wykresu,
ggplot(wykres_dane_holt1, aes(x = okres)) +
geom_line(aes(y = actual), color = "blue", linetype = "solid", size = 1) +
geom_line(aes(y = predicted), color = "red", linetype = "dashed", size = 1) +
geom_point(data = wykres_prognozowane_dane_holt1, aes(x = Time, y = Lo95), colour = "orange", alpha = 1) +
geom_point(data = wykres_prognozowane_dane_holt1, aes(x = Time, y = Lo80), colour = "orange", alpha = 0.5) +
geom_point(data = wykres_prognozowane_dane_holt1, aes(x = Time, y = Hi95), colour = "orange", alpha = 1) +
geom_point(data = wykres_prognozowane_dane_holt1, aes(x = Time, y = Hi80), colour = "orange", alpha = 0.5) +
labs(title = paste("Rzeczywiste vs przewidziane wartości metodą Holt:", var),
x = "Rok",
y = var) +
theme_minimal()
)
lista_wykresow_holt1[[nazwa_wykresu]] <- get(nazwa_wykresu)
}## [1] "Model"
## [1] "Oferty_pracy"
##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = df_ost[[var]], alpha = NULL, beta = NULL)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 1e-04
##
## Initial states:
## l = 6.0212
## b = 0.092
##
## sigma: 0.4698
##
## AIC AICc BIC
## 32.74847 37.36385 37.47066
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0046022 0.4174562 0.3327534 0.02806891 4.330741 0.9563353
## ACF1
## Training set 0.08212374
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 7.950325 7.348212 8.552438 7.029472 8.871178
## RMSPE: 5.451482
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.082 0.131 -0.177 -0.363 -0.209 -0.127 -0.082 0.285 0.131 0.172
## 11 12
## 0.072 -0.130
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.082 0.125 -0.202 -0.370 -0.142 -0.047 -0.188 0.150 0.029 -0.021 0.025
## 12
## -0.033
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -2.4244, Lag order = 2, p-value = 0.4107
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.1495, df = 1, p-value = 0.699
##
## -------------------------------------------------------------------------------
## [1] "Model"
## [1] "Prz_mie_dochód_rozp_na_1_osobę"
##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = df_ost[[var]], alpha = NULL, beta = NULL)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.1617
##
## Initial states:
## l = 6.5532
## b = 0.0719
##
## sigma: 0.0522
##
## AIC AICc BIC
## -50.73495 -46.11956 -46.01275
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.007906536 0.0463976 0.03486918 -0.1080881 0.4818398 0.5531621
## ACF1
## Training set 0.08303958
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 7.658583 7.591661 7.725504 7.556236 7.760929
## RMSPE: 0.644803
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.083 -0.348 0.122 -0.248 -0.096 0.366 -0.211 -0.397 0.129 0.132
## 11 12
## -0.007 0.075
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.083 -0.357 0.220 -0.513 0.335 -0.116 -0.189 -0.327 0.074 0.107 -0.161
## 12
## -0.148
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -1.8296, Lag order = 2, p-value = 0.6373
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.15285, df = 1, p-value = 0.6958
##
## -------------------------------------------------------------------------------
## [1] "Model"
## [1] "Inflacja"
##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = df_ost[[var]], alpha = NULL, beta = NULL)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.1486
##
## Initial states:
## l = 4.6309
## b = -5e-04
##
## sigma: 0.0251
##
## AIC AICc BIC
## -78.58045 -73.96507 -73.85826
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.005819425 0.02229743 0.0150214 0.1233809 0.3225761 0.9756551
## ACF1
## Training set 0.141507
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 4.751217 4.719057 4.783378 4.702032 4.800403
## RMSPE: 0.481398
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.142 -0.068 -0.006 0.016 0.148 0.033 -0.071 -0.113 -0.160 -0.066
## 11 12
## -0.014 -0.080
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.142 -0.090 0.017 0.009 0.149 -0.011 -0.054 -0.100 -0.149 -0.063 -0.024
## 12
## -0.066
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -0.75442, Lag order = 2, p-value = 0.9538
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.44387, df = 1, p-value = 0.5053
##
## -------------------------------------------------------------------------------
## [1] "Model"
## [1] "Wydatki_na_Oświata_i_wychowanie"
##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = df_ost[[var]], alpha = NULL, beta = NULL)
##
## Smoothing parameters:
## alpha = 0.988
## beta = 0.988
##
## Initial states:
## l = 16.5324
## b = 0.226
##
## sigma: 0.057
##
## AIC AICc BIC
## -47.41353 -42.79814 -42.69133
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.006223361 0.05063553 0.03633174 -0.0366341 0.2107245 0.5906457
## ACF1
## Training set 0.3161273
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 17.71623 17.6432 17.78927 17.60454 17.82793
## RMSPE: 0.293876
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.316 -0.199 -0.114 -0.159 -0.107 0.191 0.313 0.077 -0.175 -0.178
## 11 12
## -0.173 -0.171
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.316 -0.332 0.093 -0.265 0.056 0.161 0.176 -0.046 -0.116 -0.035 -0.128
## 12
## -0.134
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -2.8071, Lag order = 2, p-value = 0.2649
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 2.2153, df = 1, p-value = 0.1367
##
## -------------------------------------------------------------------------------rmspe_values_holt1## $Oferty_pracy
## [1] "RMSPE (holt1): 5.45148186201636"
##
## $Prz_mie_dochód_rozp_na_1_osobę
## [1] "RMSPE (holt1): 0.644803030301088"
##
## $Inflacja
## [1] "RMSPE (holt1): 0.481397965983889"
##
## $Wydatki_na_Oświata_i_wychowanie
## [1] "RMSPE (holt1): 0.293875950859098"# wyświetlenie wszystkich wykresów z listy
for (wykres in lista_wykresow_holt1) {
print(wykres)
}## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
Prognozowanie - metoda Holta (drugie podejście)
rmspe_values_holt2 <- list() #lista na wyniki RMSPE
lista_wykresow_holt2 <- list()
x_do_prognozowania_y <- list() #lista, do której trafią wyniki RMSPE
for (var in colnames(df_ost[-1])) {
# prognoza
holt_model2 <- holt(df_ost[[var]], alpha=NULL, beta=NULL, damped=TRUE)
holt_prognoza2 <- forecast(holt_model2, h = 1)
# wypisanie wynikow
print("Model")
print(var)
print(summary(holt_prognoza2))
rmspe_value <- fun_rmspe(df_ost[[var]], holt_prognoza2$fitted)
cat("RMSPE: ", rmspe_value, "\n")
rmspe_values_holt2[[var]] <- paste("RMSPE (holt2): ", rmspe_value)
if (var != "Wydatki_na_Oświata_i_wychowanie") {
x_do_prognozowania_y[[var]] <- holt_prognoza2$mean
}
cat("-------------------------------------------------------------------------------\n")
# wyliczenie reszt
print("Testy dla reszt")
res <- holt_prognoza2$resid
plot(res, ylab = paste("Reszty (", var, ")", sep = "")) #wizualizacja reszt
# Testy dla reszt
print(acf(res, type = c("correlation"), plot = TRUE, ylab = paste("ACF (", var, ")", sep = "")))
print(pacf(res, plot = TRUE, ylab = paste("PACF (", var, ")", sep = "")))
print(adf.test(res))
print(Box.test (res, lag = , type = "Ljung"))
cat("-------------------------------------------------------------------------------\n")
# utworzenie dataframe z danymi potrzebnymi do stworzenia wykresow
wykres_dane_holt2 <- data.frame(okres = c(df_log$Okres, 2023), actual = c(df_ost[[var]], NA), predicted = c(holt_prognoza2$fitted, holt_prognoza2$mean))
wykres_prognozowane_dane_holt2 <- data.frame(Time = 2023,
Forecast = as.numeric(holt_prognoza2$mean),
Lo80 = holt_prognoza2$lower[, "80%"],
Hi80 = holt_prognoza2$upper[, "80%"],
Lo95 = holt_prognoza2$lower[, "95%"],
Hi95 = holt_prognoza2$upper[, "95%"])
# narysowanie wykresow
nazwa_wykresu <- paste("wykres_holt2_", var, sep = "")
assign(nazwa_wykresu,
ggplot(wykres_dane_holt2, aes(x = okres)) +
geom_line(aes(y = actual), color = "blue", linetype = "solid", size = 1) +
geom_line(aes(y = predicted), color = "red", linetype = "dashed", size = 1) +
geom_point(data = wykres_prognozowane_dane_holt2, aes(x = Time, y = Lo95), colour = "orange", alpha = 1) +
geom_point(data = wykres_prognozowane_dane_holt2, aes(x = Time, y = Lo80), colour = "orange", alpha = 0.5) +
geom_point(data = wykres_prognozowane_dane_holt2, aes(x = Time, y = Hi95), colour = "orange", alpha = 1) +
geom_point(data = wykres_prognozowane_dane_holt2, aes(x = Time, y = Hi80), colour = "orange", alpha = 0.5) +
labs(title = paste("Rzeczywiste vs przewidziane wartości metodą Holt:", var),
x = "Rok",
y = var) +
theme_minimal()
)
lista_wykresow_holt2[[nazwa_wykresu]] <- get(nazwa_wykresu)
}## [1] "Model"
## [1] "Oferty_pracy"
##
## Forecast method: Damped Holt's method
##
## Model Information:
## Damped Holt's method
##
## Call:
## holt(y = df_ost[[var]], damped = TRUE, alpha = NULL, beta = NULL)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.0519
## phi = 0.8
##
## Initial states:
## l = 5.4545
## b = 0.9452
##
## sigma: 0.4498
##
## AIC AICc BIC
## 31.77843 38.77843 37.44507
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.06036977 0.3860708 0.3123514 -0.8830888 4.073218 0.8976999
## ACF1
## Training set 0.007248366
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 7.85665 7.280261 8.433039 6.975138 8.738161
## RMSPE: 5.041626
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.007 0.105 -0.077 -0.364 -0.222 -0.117 -0.157 0.215 0.034 0.096
## 11 12
## 0.082 -0.104
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.007 0.105 -0.079 -0.380 -0.246 -0.070 -0.212 0.035 -0.126 -0.121 -0.098
## 12
## -0.164
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -2.0053, Lag order = 2, p-value = 0.5704
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.0011646, df = 1, p-value = 0.9728
##
## -------------------------------------------------------------------------------
## [1] "Model"
## [1] "Prz_mie_dochód_rozp_na_1_osobę"
##
## Forecast method: Damped Holt's method
##
## Model Information:
## Damped Holt's method
##
## Call:
## holt(y = df_ost[[var]], damped = TRUE, alpha = NULL, beta = NULL)
##
## Smoothing parameters:
## alpha = 0.9784
## beta = 1e-04
## phi = 0.932
##
## Initial states:
## l = 6.5251
## b = 0.1064
##
## sigma: 0.0481
##
## AIC AICc BIC
## -53.14825 -46.14825 -47.48161
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0004853764 0.04131015 0.03169671 0.004345197 0.4370219 0.5028342
## ACF1
## Training set 0.01073952
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 7.635584 7.57391 7.697259 7.541261 7.729907
## RMSPE: 0.5741011
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.011 -0.411 0.141 -0.224 -0.065 0.409 -0.227 -0.402 0.195 0.137
## 11 12
## -0.064 0.065
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.011 -0.412 0.182 -0.509 0.275 -0.024 -0.215 -0.293 0.011 0.091 -0.174
## 12
## -0.150
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -1.8992, Lag order = 2, p-value = 0.6108
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.0025566, df = 1, p-value = 0.9597
##
## -------------------------------------------------------------------------------
## [1] "Model"
## [1] "Inflacja"
##
## Forecast method: Damped Holt's method
##
## Model Information:
## Damped Holt's method
##
## Call:
## holt(y = df_ost[[var]], damped = TRUE, alpha = NULL, beta = NULL)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.748
## phi = 0.893
##
## Initial states:
## l = 4.6503
## b = -0.0148
##
## sigma: 0.0249
##
## AIC AICc BIC
## -78.17513 -71.17513 -72.50850
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.005430324 0.02138107 0.01478985 0.1154753 0.3178129 0.9606154
## ACF1
## Training set 0.002627294
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 4.792961 4.76104 4.824882 4.744142 4.84178
## RMSPE: 0.461614
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.003 -0.164 0.025 -0.106 0.123 0.104 -0.044 -0.054 -0.079 -0.107
## 11 12
## -0.030 0.001
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.003 -0.164 0.027 -0.137 0.140 0.060 0.008 -0.052 -0.064 -0.128 -0.080
## 12
## -0.050
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -1.1268, Lag order = 2, p-value = 0.9019
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.00015301, df = 1, p-value = 0.9901
##
## -------------------------------------------------------------------------------
## [1] "Model"
## [1] "Wydatki_na_Oświata_i_wychowanie"
##
## Forecast method: Damped Holt's method
##
## Model Information:
## Damped Holt's method
##
## Call:
## holt(y = df_ost[[var]], damped = TRUE, alpha = NULL, beta = NULL)
##
## Smoothing parameters:
## alpha = 0.9988
## beta = 0.9988
## phi = 0.8
##
## Initial states:
## l = 16.5997
## b = 0.0914
##
## sigma: 0.0612
##
## AIC AICc BIC
## -44.03484 -37.03484 -38.36821
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01136728 0.05250638 0.04297093 0.06646406 0.2501463 0.698579
## ACF1
## Training set 0.2632309
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 17.69463 17.61624 17.77302 17.57474 17.81452
## RMSPE: 0.3047339
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.263 -0.389 -0.215 -0.088 -0.050 0.158 0.237 0.009 -0.283 -0.235
## 11 12
## -0.074 -0.051
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.263 -0.492 0.096 -0.322 0.030 0.058 0.114 -0.033 -0.199 -0.055 -0.222
## 12
## -0.175
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -3.1419, Lag order = 2, p-value = 0.1374
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 1.5359, df = 1, p-value = 0.2152
##
## -------------------------------------------------------------------------------rmspe_values_holt2## $Oferty_pracy
## [1] "RMSPE (holt2): 5.04162631094541"
##
## $Prz_mie_dochód_rozp_na_1_osobę
## [1] "RMSPE (holt2): 0.574101055892829"
##
## $Inflacja
## [1] "RMSPE (holt2): 0.461613990468139"
##
## $Wydatki_na_Oświata_i_wychowanie
## [1] "RMSPE (holt2): 0.304733891028666"# wyświetlenie wszystkich wykresów z listy
for (wykres in lista_wykresow_holt2) {
print(wykres)
}## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
Prognozowanie - ARIMA
rmspe_values_arima <- list() #lista na wyniki RMSPE
lista_wykresow_arima <- list()
for (var in colnames(df_ost[-1])) {
# prognoza
arima_model <- auto.arima(df_ost[[var]])
arima_prognoza <- forecast(arima_model, h = 1)
# wypisanie wynikow
print("Model")
print(var)
print(summary(arima_prognoza))
rmspe_value <- fun_rmspe(df_ost[[var]], arima_prognoza$fitted)
cat("RMSPE: ", rmspe_value, "\n")
rmspe_values_arima[[var]] <- paste("RMSPE (ARIMA): ", rmspe_value)
if (var == "Wydatki_na_Oświata_i_wychowanie") {
x_do_prognozowania_y[[var]] <- arima_prognoza$mean
}
cat("-------------------------------------------------------------------------------\n")
# wyliczenie reszt
print("Testy dla reszt")
res <- arima_prognoza$resid
plot(res, ylab = paste("Reszty (", var, ")", sep = "")) #wizualizacja reszt
# Testy dla reszt
print(acf(res, type = c("correlation"), plot = TRUE, ylab = paste("ACF (", var, ")", sep = "")))
print(pacf(res, plot = TRUE, ylab = paste("PACF (", var, ")", sep = "")))
print(adf.test(res))
print(Box.test (res, lag = , type = "Ljung"))
cat("-------------------------------------------------------------------------------\n")
# utworzenie dataframe z danymi potrzebnymi do stworzenia wykresow
wykres_dane_arima <- data.frame(okres = c(df_log$Okres, 2023), actual = c(df_ost[[var]], NA), predicted = c(arima_prognoza$fitted, arima_prognoza$mean))
wykres_prognozowane_dane_arima <- data.frame(Time = 2023,
Forecast = as.numeric(arima_prognoza$mean),
Lo80 = arima_prognoza$lower[, "80%"],
Hi80 = arima_prognoza$upper[, "80%"],
Lo95 = arima_prognoza$lower[, "95%"],
Hi95 = arima_prognoza$upper[, "95%"])
# narysowanie wykresow
nazwa_wykresu <- paste("wykres_arima_", var, sep = "")
assign(nazwa_wykresu,
ggplot(wykres_dane_arima, aes(x = okres)) +
geom_line(aes(y = actual), color = "blue", linetype = "solid", size = 1) +
geom_line(aes(y = predicted), color = "red", linetype = "dashed", size = 1) +
geom_point(data = wykres_prognozowane_dane_arima, aes(x = Time, y = Lo95), colour = "orange", alpha = 1) +
geom_point(data = wykres_prognozowane_dane_arima, aes(x = Time, y = Lo80), colour = "orange", alpha = 0.5) +
geom_point(data = wykres_prognozowane_dane_arima, aes(x = Time, y = Hi95), colour = "orange", alpha = 1) +
geom_point(data = wykres_prognozowane_dane_arima, aes(x = Time, y = Hi80), colour = "orange", alpha = 0.5) +
labs(title = paste("Rzeczywiste vs przewidziane wartości metodą Holt:", var),
x = "Rok",
y = var) +
theme_minimal()
)
lista_wykresow_arima[[nazwa_wykresu]] <- get(nazwa_wykresu)
}## [1] "Model"
## [1] "Oferty_pracy"
##
## Forecast method: ARIMA(0,1,0)
##
## Model Information:
## Series: df_ost[[var]]
## ARIMA(0,1,0)
##
## sigma^2 = 0.1933: log likelihood = -10.75
## AIC=23.49 AICc=23.74 BIC=24.38
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.09190921 0.4278829 0.3299554 1.159898 4.285444 0.9482939
## ACF1
## Training set 0.05922519
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 7.858254 7.294874 8.421634 6.996639 8.71987
## RMSPE: 5.587642
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.059 0.124 -0.201 -0.343 -0.198 -0.125 -0.069 0.288 0.123 0.165
## 11 12
## 0.058 -0.130
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.059 0.121 -0.218 -0.354 -0.143 -0.080 -0.199 0.150 0.021 -0.033 0.036
## 12
## -0.028
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -2.4821, Lag order = 2, p-value = 0.3887
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.077752, df = 1, p-value = 0.7804
##
## -------------------------------------------------------------------------------
## [1] "Model"
## [1] "Prz_mie_dochód_rozp_na_1_osobę"
##
## Forecast method: ARIMA(0,1,0) with drift
##
## Model Information:
## Series: df_ost[[var]]
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## 0.0548
## s.e. 0.0108
##
## sigma^2 = 0.002239: log likelihood = 29.9
## AIC=-55.8 AICc=-55 BIC=-54.02
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0003457726 0.04475506 0.03182293 0.01176775 0.4409394 0.5048366
## ACF1
## Training set 0.2046413
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 7.665834 7.605198 7.72647 7.573099 7.758568
## RMSPE: 0.6219761
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.205 -0.217 0.137 -0.197 -0.072 0.324 -0.164 -0.354 0.068 0.091
## 11 12
## -0.029 0.011
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.205 -0.270 0.280 -0.457 0.425 -0.152 -0.114 -0.274 0.217 0.076 -0.167
## 12
## -0.157
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -1.8284, Lag order = 2, p-value = 0.6378
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.9283, df = 1, p-value = 0.3353
##
## -------------------------------------------------------------------------------
## [1] "Model"
## [1] "Inflacja"
##
## Forecast method: ARIMA(1,0,0) with non-zero mean
##
## Model Information:
## Series: df_ost[[var]]
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.8224 4.6496
## s.e. 0.1804 0.0283
##
## sigma^2 = 0.000579: log likelihood = 44.35
## AIC=-82.7 AICc=-81.1 BIC=-79.86
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0008527541 0.02276071 0.01523295 0.01572059 0.3272061 0.9893951
## ACF1
## Training set 0.2090452
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 4.720094 4.689257 4.750931 4.672932 4.767255
## RMSPE: 0.4914001
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.209 0.024 -0.011 -0.003 0.109 -0.059 -0.159 -0.147 -0.124 0.005
## 11 12
## 0.055 -0.046##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.209 -0.020 -0.012 0.002 0.115 -0.112 -0.133 -0.088 -0.080 0.029 0.067
## 12
## -0.050## Warning in adf.test(res): p-value greater than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -0.082499, Lag order = 2, p-value = 0.99
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.96868, df = 1, p-value = 0.325
##
## -------------------------------------------------------------------------------
## [1] "Model"
## [1] "Wydatki_na_Oświata_i_wychowanie"
##
## Forecast method: ARIMA(2,1,0) with drift
##
## Model Information:
## Series: df_ost[[var]]
## ARIMA(2,1,0) with drift
##
## Coefficients:
## ar1 ar2 drift
## 1.1245 -0.5609 0.0556
## s.e. 0.2314 0.2372 0.0273
##
## sigma^2 = 0.002333: log likelihood = 29.97
## AIC=-51.95 AICc=-48.87 BIC=-48.38
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001708513 0.04291899 0.03278296 -0.009613184 0.190417 0.532953
## ACF1
## Training set -0.1066372
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 17.70815 17.64625 17.77006 17.61348 17.80283
## RMSPE: 0.2490911
## -------------------------------------------------------------------------------
## [1] "Testy dla reszt"
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 -0.107 -0.273 0.354 0.013 -0.152 -0.022 0.081 -0.161 -0.320 0.025
## 11 12
## 0.044 -0.246
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## -0.107 -0.287 0.315 -0.005 0.025 -0.169 0.050 -0.198 -0.313 -0.193 -0.048
## 12
## -0.138
##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -1.5801, Lag order = 2, p-value = 0.7323
## alternative hypothesis: stationary
##
##
## Box-Ljung test
##
## data: res
## X-squared = 0.25207, df = 1, p-value = 0.6156
##
## -------------------------------------------------------------------------------rmspe_values_arima## $Oferty_pracy
## [1] "RMSPE (ARIMA): 5.58764239605644"
##
## $Prz_mie_dochód_rozp_na_1_osobę
## [1] "RMSPE (ARIMA): 0.62197609729362"
##
## $Inflacja
## [1] "RMSPE (ARIMA): 0.491400118884011"
##
## $Wydatki_na_Oświata_i_wychowanie
## [1] "RMSPE (ARIMA): 0.249091052680956"# wyświetlenie wszystkich wykresów z listy
for (wykres in lista_wykresow_arima) {
print(wykres)
}## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 1 row containing missing values (`geom_line()`).
Prognozowanie zmiennej objaśniającej
Wydobycie wyników RMSPE z każdej prognozy zmiennych objaśniających
# funkcja, która wydobywa wyniki RMSPE
extract_rmspe <- function(text) {
gsub(".*: ([0-9.]+).*", "\\1", text)
}
# wydobycie wyników RMSPE
RMSPE_wyniki_razem <- as.data.frame(cbind(rmspe_values_lm, rmspe_values_lm2, rmspe_values_ma, rmspe_values_holt1, rmspe_values_holt2, rmspe_values_arima))
RMSPE_wyniki_razem <- as.data.frame(lapply(RMSPE_wyniki_razem, extract_rmspe))
RMSPE_wyniki_razem## rmspe_values_lm rmspe_values_lm2
## 1 RMSPE (lm): 5.8143692677273 RMSPE (lm): 5.60037960941203
## 2 RMSPE (lm): 0.861617734292072 RMSPE (lm): 0.610283665221309
## 3 RMSPE (lm): 0.581073268277916 RMSPE (lm): 0.470107826683551
## 4 RMSPE (lm): 0.539973489648822 RMSPE (lm): 0.538356643913684
## rmspe_values_ma rmspe_values_holt1
## 1 RMSPE (MA): 6.90563841214671 RMSPE (holt1): 5.45148186201636
## 2 RMSPE (MA): 0.95740707253437 RMSPE (holt1): 0.644803030301088
## 3 RMSPE (MA): 0.613652739293226 RMSPE (holt1): 0.481397965983889
## 4 RMSPE (MA): 0.625956299173227 RMSPE (holt1): 0.293875950859098
## rmspe_values_holt2 rmspe_values_arima
## 1 RMSPE (holt2): 5.04162631094541 RMSPE (ARIMA): 5.58764239605644
## 2 RMSPE (holt2): 0.574101055892829 RMSPE (ARIMA): 0.62197609729362
## 3 RMSPE (holt2): 0.461613990468139 RMSPE (ARIMA): 0.491400118884011
## 4 RMSPE (holt2): 0.304733891028666 RMSPE (ARIMA): 0.249091052680956Na podstawie otrzymanych wartości błędu RMSPE wybraliśmy prognozowane wartości zmiennych objaśniających do finalnej prognozy zmiennej objaśnianej czyli liczby zarejestrowanych osób bezrobotnych na Pomorzu w roku 2023.
Dla zmiennej Oferty Pracy najmniejszy błąd RMSPE wynoszący 5,04% miała druga prognoza metodą Holta.
Dla zmiennej Przeciętny Miesięczny Dochód Rozporządzalny na 1 osobę najmniejszy błąd RMSPE również miała druga prognoza metodą Holta,wartość tego błędu wyniosła 0,57%
Dla zmiennej Inflacja również najlepsza okazała się druga prognoza metodą Holta, w tym przypadku błąd RMSPE wyniósł 0,46%
Dla zmiennej Wydatki na Oświatę i Wychowanie najmniejszą wartość błędu RMSPE miała prognoza metodą ARIMA. Błąd ten wyniósł 0,25%
Utworzenie zbioru danych przeznaczonego do zaprognozowania zmiennej objaśnianej
# utworzenie zbioru danych przeznaczonego do zaprognozowania zmiennej objaśnianej
df_y <- subset(df_log, select = c(Okres, Bezrobotni_zarejestrowani, Oferty_pracy, Prz_mie_dochód_rozp_na_1_osobę, Inflacja, Wydatki_na_Oświata_i_wychowanie))
df_y # zbiór danych przeznaczony do zaprognozowania zmiennej objaśnianej## # A tibble: 19 × 6
## Okres Bezrobotni_zarejestrowani Oferty_pracy Prz_mie_dochód_rozp_n…¹ Inflacja
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2004 12.1 6.12 6.62 4.64
## 2 2005 12.0 7.06 6.69 4.62
## 3 2006 11.7 7.37 6.81 4.61
## 4 2007 11.4 8.30 6.94 4.63
## 5 2008 11.1 7.67 7.01 4.64
## 6 2009 11.5 7.34 7.06 4.63
## 7 2010 11.6 7.41 7.13 4.63
## 8 2011 11.6 7.04 7.16 4.65
## 9 2012 11.6 7.02 7.21 4.64
## 10 2013 11.6 7.41 7.27 4.61
## 11 2014 11.5 7.71 7.23 4.60
## 12 2015 11.3 8.30 7.23 4.60
## 13 2016 11.1 8.40 7.35 4.60
## 14 2017 10.8 8.34 7.41 4.63
## 15 2018 10.7 8.15 7.46 4.62
## 16 2019 10.6 7.92 7.53 4.62
## 17 2020 10.9 7.84 7.50 4.64
## 18 2021 10.8 8.21 7.52 4.65
## 19 2022 10.7 7.86 7.61 4.74
## # ℹ abbreviated name: ¹Prz_mie_dochód_rozp_na_1_osobę
## # ℹ 1 more variable: Wydatki_na_Oświata_i_wychowanie <dbl>#x_do_prognozowania_y <- as.data.frame(x_do_prognozowania_y)
# złączenie zbioru danych przeznaczonego do zaprognozowania zmiennej objaśnianej z zaprognozowanymi wartościami zmiennych objaśniających na 2023 rok
df_prognoza_y <- rbind(df_y, c(Okres=2023, Bezrobotni_zarejestrowani=NA, Oferty_pracy=x_do_prognozowania_y$Oferty_pracy, Prz_mie_dochód_rozp_na_1_osobę=x_do_prognozowania_y$Prz_mie_dochód_rozp_na_1_osobę, Inflacja=x_do_prognozowania_y$Inflacja, Wydatki_na_Oświata_i_wychowanie=x_do_prognozowania_y$Wydatki_na_Oświata_i_wychowanie))
df_prognoza_y # zaaktualizowany zbiór danych przeznaczony do zaprognozowania zmiennej objaśnianej## # A tibble: 20 × 6
## Okres Bezrobotni_zarejestrowani Oferty_pracy Prz_mie_dochód_rozp_n…¹ Inflacja
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2004 12.1 6.12 6.62 4.64
## 2 2005 12.0 7.06 6.69 4.62
## 3 2006 11.7 7.37 6.81 4.61
## 4 2007 11.4 8.30 6.94 4.63
## 5 2008 11.1 7.67 7.01 4.64
## 6 2009 11.5 7.34 7.06 4.63
## 7 2010 11.6 7.41 7.13 4.63
## 8 2011 11.6 7.04 7.16 4.65
## 9 2012 11.6 7.02 7.21 4.64
## 10 2013 11.6 7.41 7.27 4.61
## 11 2014 11.5 7.71 7.23 4.60
## 12 2015 11.3 8.30 7.23 4.60
## 13 2016 11.1 8.40 7.35 4.60
## 14 2017 10.8 8.34 7.41 4.63
## 15 2018 10.7 8.15 7.46 4.62
## 16 2019 10.6 7.92 7.53 4.62
## 17 2020 10.9 7.84 7.50 4.64
## 18 2021 10.8 8.21 7.52 4.65
## 19 2022 10.7 7.86 7.61 4.74
## 20 2023 NA 7.86 7.64 4.79
## # ℹ abbreviated name: ¹Prz_mie_dochód_rozp_na_1_osobę
## # ℹ 1 more variable: Wydatki_na_Oświata_i_wychowanie <dbl>Prognozowanie zmiennej objaśnianej (Bezrobocie zarejestrowane) - metoda Holta
rmspe_values_y <- list() #lista na wyniki RMSPE
lista_wykresow_y <- list()
# prognoza
y_model <- holt(df_prognoza_y[["Bezrobotni_zarejestrowani"]], alpha=NULL, beta=NULL, damped=TRUE)## Warning in ets(x, "AAN", alpha = alpha, beta = beta, phi = phi, damped =
## damped, : Missing values encountered. Using longest contiguous portion of time
## seriesy_prognoza <- forecast(y_model, h = 1)
# wypisanie wynikow
print("Model")## [1] "Model"print("Bezrobotni_zarejestrowani")## [1] "Bezrobotni_zarejestrowani"print(summary(y_prognoza))##
## Forecast method: Damped Holt's method
##
## Model Information:
## Damped Holt's method
##
## Call:
## holt(y = df_prognoza_y[["Bezrobotni_zarejestrowani"]], damped = TRUE,
##
## Call:
## alpha = NULL, beta = NULL)
##
## Smoothing parameters:
## alpha = 0.9992
## beta = 1e-04
## phi = 0.9417
##
## Initial states:
## l = 12.2188
## b = -0.1372
##
## sigma: 0.212
##
## AIC AICc BIC
## 3.201577 10.201577 8.868210
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.002951117 0.1819985 0.1405446 -0.04068567 1.248358 0.8194821
## ACF1
## Training set 0.1822039
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 20 10.61372 10.34201 10.88544 10.19817 11.02928rmspe_value <- fun_rmspe(df_prognoza_y[["Bezrobotni_zarejestrowani"]][-20], y_prognoza$fitted)
cat("RMSPE: ", rmspe_value, "\n")## RMSPE: 1.61111rmspe_values_y[["Bezrobotni_zarejestrowani"]] <- paste("RMSPE: ", rmspe_value)
cat("-------------------------------------------------------------------------------\n")## -------------------------------------------------------------------------------# wyliczenie reszt
print("Testy dla reszt")## [1] "Testy dla reszt"res <- y_prognoza$resid
plot(res, ylab = paste("Reszty (Bezrobotni_zarejestrowani)", sep = "")) #wizualizacja reszt
# Testy dla reszt
print(acf(res, type = c("correlation"), plot = TRUE, ylab = paste("ACF (Bezrobotni_zarejestrowani)", sep = "")))
##
## Autocorrelations of series 'res', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.182 -0.083 -0.172 -0.164 -0.298 -0.231 -0.016 0.034 0.126 0.111
## 11 12
## 0.283 -0.169print(pacf(res, plot = TRUE, ylab = paste("PACF (Bezrobotni_zarejestrowani)", sep = "")))
##
## Partial autocorrelations of series 'res', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.182 -0.120 -0.140 -0.122 -0.299 -0.227 -0.095 -0.169 -0.060 -0.111 0.150
## 12
## -0.366print(adf.test(res))##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -2.3194, Lag order = 2, p-value = 0.4507
## alternative hypothesis: stationaryprint(Box.test (res, lag = , type = "Ljung"))##
## Box-Ljung test
##
## data: res
## X-squared = 0.73589, df = 1, p-value = 0.391cat("-------------------------------------------------------------------------------\n")## -------------------------------------------------------------------------------# utworzenie dataframe z danymi potrzebnymi do stworzenia wykresu
wykres_dane_y <- data.frame(okres = c(df_prognoza_y$Okres[-20], 2023), actual = c(df_prognoza_y[["Bezrobotni_zarejestrowani"]][-20], NA), predicted = c(y_prognoza$fitted, y_prognoza$mean))
wykres_prognozowane_dane_y <- data.frame(Time = 2023,
Forecast = as.numeric(y_prognoza$mean),
Lo80 = y_prognoza$lower[, "80%"],
Hi80 = y_prognoza$upper[, "80%"],
Lo95 = y_prognoza$lower[, "95%"],
Hi95 = y_prognoza$upper[, "95%"])
# narysowanie wykresu
wykres_y <- ggplot(wykres_dane_y, aes(x = okres)) +
geom_line(aes(y = actual), color = "blue", linetype = "solid", size = 1) +
geom_line(aes(y = predicted), color = "red", linetype = "dashed", size = 1) +
geom_point(data = wykres_prognozowane_dane_y, aes(x = Time, y = Lo95), colour = "orange", alpha = 1) +
geom_point(data = wykres_prognozowane_dane_y, aes(x = Time, y = Lo80), colour = "orange", alpha = 0.5) +
geom_point(data = wykres_prognozowane_dane_y, aes(x = Time, y = Hi95), colour = "orange", alpha = 1) +
geom_point(data = wykres_prognozowane_dane_y, aes(x = Time, y = Hi80), colour = "orange", alpha = 0.5) +
labs(title = paste("Rzeczywiste vs przewidziane wartości metodą Holt:", "Bezrobotni_zarejestrowani"),
x = "Rok",
y = "Bezrobotni_zarejestrowani") +
theme_minimal()
wykres_y## Warning: Removed 1 row containing missing values (`geom_line()`).
rmspe_values_y## $Bezrobotni_zarejestrowani
## [1] "RMSPE: 1.61110971588003"#wypisanie wyników
expv1<-exp(y_prognoza$mean[1])
cat("Wynik prognozy metodą Holta:", y_prognoza$mean[1],"\n")## Wynik prognozy metodą Holta: 10.61372cat("Wynik prognozy metodą Holta po odlogarytmowaniu:", expv1)## Wynik prognozy metodą Holta po odlogarytmowaniu: 40689.41Średnia oczekiwana wartość otrzymana na podstawie prognozy wynosi 10.61372. Jest to wartość zlogarytmowana, co oznacza, że po odlogarytmowaniu wartość wyniesie około 40689. Oznacza to, że przewidywana wartość liczby osób bezrobotnych w 2023 roku na Pomorzu wyniesie około 40689 osób - wartość tą otrzymano przy użyciu funkcji exp umożliwiającej odlogarytmowanie wartości. Pozostałe przedstawione w powyższej tabeli wartości logarytmiczne oznaczają odpowiednio: Low80- 20% obserwacji jest równych tej wartości lub mniejszych, Hi80- 20% obserwacji jest równych tej wartości lub większych, Low95- 5% obserwacji jest rónych tej wartości lub mniejszych, Hi95- 5% obserwacji jest równych tej wartości lub większych.
Prognoza wartości zmiennej objaśnianej w 2023 roku na podstawie modelu
# utworzenie dataframe z wartościami x do prognozy y
prognozaost <- data.frame(Oferty_pracy=x_do_prognozowania_y$Oferty_pracy, Prz_mie_dochód_rozp_na_1_osobę=x_do_prognozowania_y$Prz_mie_dochód_rozp_na_1_osobę, Inflacja=x_do_prognozowania_y$Inflacja, Wydatki_na_Oświata_i_wychowanie=x_do_prognozowania_y$Wydatki_na_Oświata_i_wychowanie)
# prognoza y
wartosc_y <- predict(model_log, newdata = prognozaost)
cat("Wynik prognozy na podstawie modelu:", wartosc_y,"\n") #wynik: 10,31524## Wynik prognozy na podstawie modelu: 10.31524# odlogarytmowanie wyniku
expv2 <- exp(wartosc_y)
cat("Wynik prognozy na podstawie modelu po odlogarytmowaniu:", expv2) #wynik: 30189.22## Wynik prognozy na podstawie modelu po odlogarytmowaniu: 30189.22# utworzenie dataframe z danymi potrzebnymi do stworzenia wykresu
wykres_dane_model <- data.frame(okres = c(df_prognoza_y$Okres[-20], 2023), actual = c(df_prognoza_y[["Bezrobotni_zarejestrowani"]][-20], NA), predicted = c(fitted(model_log), wartosc_y))
wykres_y <- ggplot(wykres_dane_model, aes(x = okres)) +
geom_line(aes(y = actual), color = "blue", linetype = "solid", size = 1) +
geom_line(aes(y = predicted), color = "red", linetype = "dashed", size = 1)+
labs(title = paste("Rzeczywiste vs przewidziane wartości modelu:", "Bezrobotni_zarejestrowani"),
x = "Rok",
y = "Bezrobotni_zarejestrowani")+
theme_minimal()
wykres_y## Warning: Removed 1 row containing missing values (`geom_line()`).
W przypadku prognozy na podstawie modelu przewidywana liczba zarejestrowanych osób bezrobotnych na Pomorzu wynosi 30 189.
Podsumowanie
Zaprognozowanie ilości zarejestrowanych osób bezrobotnych na Pomorzu w roku 2023 wymagało na początku wybrania zmiennych które potencjalnie mogły wpływać na kształtowanie się tej zmiennej. Ostatecznie nasza grupa wybrała w tym celu 8 różnych zmiennych.
Podczas tworzenia modelu ostatecznie statystycznie istotne okazały się 4 zmienne: Oferty pracy na Pomorzu, przeciętny miesięczny dochód rozporządzalny na 1 osobę na Pomorzu, odnotowany poziom inflacji na Pomorzu oraz wydatki na oświatę i wychowanie na Pomorzu. Następnie należało wykonać testy statystyczne w celu sprawdzenia poprawności zbudowanego modelu.
Kolejnym etapem pracy były prognozy zmiennych objaśniających na rok 2023 potrzebne do zaprognozowania liczby zarejestrowanych osób bezrobotnych w tym roku. W celu uzyskania jak najlepszych wyników każdą z czterech zmiennych objaśniających zaprognozowaliśmy za pomocą 6 różnych metod.
Kolejnym krokiem było porównanie błędów RMSPE dla każdej zaprognozowanej zmiennej objaśniającej. Na tej podstawie wybrane zostały wartości tych zmiennych do prognozy zmiennej objaśnianej (Bezrobotni_zarejestrowani). Do owej prognozy wpierw zastosowaliśmy metodę Holta. Wynikło z niej, że w 2023 roku na Pomorzu osób bezrobotnych było około 40689.
Kolejno zaprognozowaliśmy tę zmienną na podstawie modelu utworzonego na początku raportu. Z tej prognozy wynikło natomiast, że w 2023 roku na Pomorzu osób bezrobotnych było 30 189.