Skup mleka
TP
2020-12-29
Przedmiotem prognozy jest miesięczna wielkość skupu mleka. Źródłem danych do zbudowania prognozy jest Biuletyn Statystyczny GUS 10/2020; https://stat.gov.pl/obszary-tematyczne/inne-opracowania/ )
Wykorzystywanym narzędziem obliczeniowym są pakiety funkcji forecast(w wersji 8.13 lub wyższej) oraz fpp2:
library ("forecast")
library ("fpp2")## ── Attaching packages ────────────────────────────────────────────── fpp2 2.4 ──## ✔ ggplot2 3.3.2 ✔ expsmooth 2.3
## ✔ fma 2.4## ## pokaż wersje pakietów
sessionInfo()## R version 3.5.3 (2019-03-11)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] expsmooth_2.3 fma_2.4 ggplot2_3.3.2 fpp2_2.4 forecast_8.13
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.3 urca_1.3-0 pillar_1.4.6 compiler_3.5.3
## [5] tseries_0.10-47 tools_3.5.3 xts_0.12-0 digest_0.6.23
## [9] nlme_3.1-151 evaluate_0.14 lifecycle_0.2.0 tibble_3.0.3
## [13] gtable_0.3.0 lattice_0.20-41 pkgconfig_2.0.3 rlang_0.4.9
## [17] cli_1.1.0 rstudioapi_0.13 curl_4.3 yaml_2.2.0
## [21] parallel_3.5.3 xfun_0.19 withr_2.1.2 stringr_1.4.0
## [25] knitr_1.26 vctrs_0.3.4 lmtest_0.9-38 grid_3.5.3
## [29] nnet_7.3-14 glue_1.4.2 R6_2.4.1 rmarkdown_2.6
## [33] purrr_0.3.3 TTR_0.24.2 magrittr_1.5 scales_1.1.0
## [37] ellipsis_0.3.0 htmltools_0.4.0 assertthat_0.2.1 quantmod_0.4.17
## [41] timeDate_3043.102 colorspace_1.4-1 fracdiff_1.5-1 quadprog_1.5-8
## [45] stringi_1.4.3 munsell_0.5.0 crayon_1.3.4 zoo_1.8-7Plik mleko.csv zawiera miesięczne dane nt skupu mleka w Polsce w okresie styczeń 2010–październik 2020 (tabela 47: skup ważniejszych produktów rolnych; https://stat.gov.pl/obszary-tematyczne/inne-opracowania/)
## załadowanie danych z pliku CSV zamiana na typ TS (szereg czasowy)
t <- read.csv("mleko.csv", dec=",", sep = ';', header=T, na.string="NA");
t <-ts(t, start=c(2010, 1), frequency=12)
## sprawdzenie ostatniej obserwacji (powinna być 2020/10)
end(t)## [1] 2020 10Wizualne rozpoznanie natury prognozowanego zjawiska (czy dane są kompletne, czy występują obserwacje nietypowe. Czy występuje trend i sezonowść.)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Analizując wykres dostrzegamy wyraźne wahania sezonowe oraz trend. Nie widać żadnych braków danych, obserwacji nietypowych czy załamań trendu. Kolejne dwa wykresy pozwolają na bardziej dokładną ocenę sezonowości i trendu:
Wielkość skup mleka spada zimą a jest najwyższa w miesiącach maj–lipiec. Ostatni wykres pokazuje wyraźnie, że sezonowość skupu ma charakter addytywny (stała amplituda wahań) a trend jest zbliżony do liniowego
W analizowanym okresie średni skup mleka wynosił 871.4907692 (wartość minimalna 625; wartość maksymalna 1074.9)
Założenia prognostyczne
Prognoza zostanie zbudowana przy założeniu że zamawiający zażądał prognozy na 6 miesięcy. W związku z tym długość zbioru testowego ustalamy na okres 6 miesięcy.
## Koniec zbioru uczącego 2020 4## Początek zbioru testowego 2020 5Prognozowanie za pomocą trendu liniowego
Oszacowanie modelu
##
## Call:
## tslm(formula = tl ~ trend + season)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.393 -13.805 1.808 14.641 36.312
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 694.39498 6.77631 102.474 < 2e-16 ***
## trend 2.49025 0.05034 49.470 < 2e-16 ***
## season2 -59.24479 8.54291 -6.935 2.84e-10 ***
## season3 27.79223 8.54335 3.253 0.001513 **
## season4 29.36562 8.54409 3.437 0.000829 ***
## season5 94.92049 8.75430 10.843 < 2e-16 ***
## season6 61.21025 8.75387 6.992 2.14e-10 ***
## season7 71.58000 8.75372 8.177 5.25e-13 ***
## season8 49.64975 8.75387 5.672 1.14e-07 ***
## season9 -1.34049 8.75430 -0.153 0.878579
## season10 -11.98074 8.75503 -1.368 0.173937
## season11 -63.23098 8.75604 -7.221 6.84e-11 ***
## season12 -20.60123 8.75734 -2.352 0.020412 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.03 on 111 degrees of freedom
## Multiple R-squared: 0.965, Adjusted R-squared: 0.9612
## F-statistic: 254.8 on 12 and 111 DF, p-value: < 2.2e-16## ME RMSE MAE MPE MAPE MASE
## Training set 2.751006e-15 18.95528 15.72372 -0.05654159 1.844314 0.4956117
## ACF1
## Training set 0.7472995Ocena reszt
##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals from Linear regression model
## LM test = 86.627, df = 24, p-value = 5.174e-09
##
## Box-Ljung test
##
## data: res_lm
## X-squared = 70.938, df = 1, p-value < 2.2e-16Wnioski: Model jest dobrze dopasowany do danych (R > 90?) ale występuje autokorelacja składnika losowego
Wyznacznie prognoz na 6 miesięcy i porównanie z wartościami zbioru testowgo
## ME RMSE MAE MPE MAPE MASE
## Training set 2.751006e-15 18.95528 15.72372 -0.05654159 1.844314 0.4956117
## Test set -3.752458e+01 38.25219 37.52458 -3.71222640 3.712226 1.1827747
## ACF1 Theil's U
## Training set 0.7472995 NA
## Test set 0.3251125 1.252958Prognozowanie za pomocą wygładzania wykładniczego
Oszacowanie modelu
## ETS(M,A,A)
##
## Call:
## ets(y = tl)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 1e-04
## gamma = 1e-04
##
## Initial states:
## l = 704.7449
## b = 2.8432
## s = -37.078 -78.6789 -26.7153 -14.825 35.1847 55.9225
## 46.1088 81.7877 17.0079 13.942 -76.0851 -16.5713
##
## sigma: 0.0166
##
## AIC AICc BIC
## 1273.366 1279.140 1321.311
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.4209157 13.46273 10.77963 -0.06212104 1.241475 0.3397738
## ACF1
## Training set 0.01570334## ME RMSE MAE MPE MAPE MASE
## Training set -0.4209157 13.46273 10.77963 -0.06212104 1.241475 0.3397738
## ACF1
## Training set 0.01570334Ocena reszt
##
## Ljung-Box test
##
## data: Residuals from ETS(M,A,A)
## Q* = 40.338, df = 8, p-value = 2.771e-06
##
## Model df: 16. Total lags used: 24
##
## Box-Ljung test
##
## data: res_es
## X-squared = 0.0068184, df = 1, p-value = 0.9342Wyznacznie prognoz na 6 miesięcy i porównanie z wartościami zbioru testowgo
## ME RMSE MAE MPE MAPE MASE
## Training set -0.4209157 13.46273 10.77963 -0.06212104 1.241475 0.3397738
## Test set -26.0203135 27.19583 26.02031 -2.58318128 2.583181 0.8201602
## ACF1 Theil's U
## Training set 0.01570334 NA
## Test set 0.36428216 0.9051521Wnioski:
Arima
## Series: tl
## ARIMA(2,0,3)(0,1,2)[12] with drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 sma1 sma2 drift
## 1.4719 -0.7928 -0.8004 0.3435 0.3039 -0.8846 0.2945 2.4476
## s.e. 0.1176 0.0970 0.1367 0.1151 0.1066 0.1237 0.1250 0.1131
##
## sigma^2 estimated as 160.3: log likelihood=-444.28
## AIC=906.55 AICc=908.32 BIC=931.02
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04190468 11.59625 8.700261 0.005871553 0.9936613 0.2742322
## ACF1
## Training set 0.02026481## ME RMSE MAE MPE MAPE MASE
## Training set 0.04190468 11.59625 8.700261 0.005871553 0.9936613 0.2742322
## ACF1
## Training set 0.02026481Ocena reszt
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,3)(0,1,2)[12] with drift
## Q* = 18.983, df = 16, p-value = 0.2696
##
## Model df: 8. Total lags used: 24
##
## Box-Ljung test
##
## data: res_aa
## X-squared = 0.052164, df = 1, p-value = 0.8193Wyznacznie prognoz na 6 miesięcy i porównanie z wartościami zbioru testowgo
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04190468 11.59625 8.700261 0.005871553 0.9936613 0.2742322
## Test set -17.10141313 18.63020 17.101413 -1.705137900 1.7051379 0.5390365
## ACF1 Theil's U
## Training set 0.02026481 NA
## Test set 0.48214641 0.6121744Wnioski
Porównanie prognoz
Ponieważ szereg wykazuje się trendem i sezonowścią porównamy wyniki do naiwnych modeli snaive (sezonowość) oraz rwf (trend)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 28.45089 38.43489 31.72589 3.256443 3.642769 1.0000000 0.7765559
## Test set 22.96667 25.03344 22.96667 2.252836 2.252836 0.7239092 -0.3385134
## Theil's U
## Training set NA
## Test set 0.8052082## ME RMSE MAE MPE MAPE MASE
## Training set 2.778049 47.74255 39.85935 0.1705065 4.669591 1.256367
## Test set -3.516667 39.99039 33.91667 -0.5002244 3.365275 1.069053
## ACF1 Theil's U
## Training set -0.2989360 NA
## Test set 0.4266543 1.176029Zestawienie ocen dopasowania (uporządkowanych wg wielkości RMSE)
## ME RMSE MAE MPE MAPE MASE
## arima 4.190468e-02 11.59625 8.700261 0.005871553 0.9936613 0.2742322
## es -4.209157e-01 13.46273 10.779628 -0.062121038 1.2414748 0.3397738
## arima/t -1.710141e+01 18.63020 17.101413 -1.705137900 1.7051379 0.5390365
## lm 2.751006e-15 18.95528 15.723723 -0.056541593 1.8443136 0.4956117
## sna/t 2.296667e+01 25.03344 22.966667 2.252836197 2.2528362 0.7239092
## es/t -2.602031e+01 27.19583 26.020314 -2.583181278 2.5831813 0.8201602
## lm/t -3.752458e+01 38.25219 37.524583 -3.712226403 3.7122264 1.1827747
## sna 2.845089e+01 38.43489 31.725893 3.256443296 3.6427687 1.0000000
## ref/t -3.516667e+00 39.99039 33.916667 -0.500224429 3.3652746 1.0690532
## rwf 2.778049e+00 47.74255 39.859350 0.170506483 4.6695910 1.2563665
## ACF1 Theil's U
## arima 0.02026481 NA
## es 0.01570334 NA
## arima/t 0.48214641 0.6121744
## lm 0.74729955 NA
## sna/t -0.33851337 0.8052082
## es/t 0.36428216 0.9051521
## lm/t 0.32511250 1.2529582
## sna 0.77655588 NA
## ref/t 0.42665435 1.1760293
## rwf -0.29893598 NAPrzyjmując jako kryterium najniższą wartość RMSE/MAPE w zbiorze testowym wybiermy model ARIMA jako najlepiej prognozujący skup mleka
Porównanie wyników prognozowania na wykresie:
Bardziej szczegółowy wykres
Wszystkie prognozy są systematycznie zawyżone się okazuje. Być może jest to spowodowane nadzwyczajnym spadkiem skupu związanym z epidemię COVID19.