Prognozowanie ilości przestępstw w gospodarce odpadami

Rządowy projekt ustawy o zmianie niektórych ustaw w celu przeciwdziałania przestępczości środowiskowej znajduje się na ostatniej prostej. 22 lipca br. ustawę przekazano Prezydentowi i Marszałkowi Senatu. Należy zatem oczekiwać, że wkrótce zacznie ona funkcjonować w polskim systemie prawnym.

Podstawowym celem projektowanej regulacji jest ochrona środowiska naturalnego przed szkodliwymi działaniami, spełniającymi znamiona przestępstw lub wykroczeń przeciwko środowisku. Listę proponowanych zmian w systemie prawnym przedstawiono na łamach portalu Sozosfera, stąd nie będę ich w tym miejscu powtarzał. Zachęcam do przeczytania artykułu na portalu Sozosfery.

Przestępstwa przeciwko środowisku, popełnione z artykułów 181-188 kodeksu karnego, policja odnotowuje w swoich statystykach. Te upublicznione dotychczas kończą się na roku 2020. Mimo wysłanej prośby do prowadzących portal o aktualizację danych za rok 2021 jeszcze aktualizacji się nie doczekałem. Zatem daję minus ;). Ale chcę też pochwalić. Na plus zaliczę wdrożenie Krajowej Mapy Zagrożeń Bezpieczeństwa, która przyspiesza i ułatwia komunikację ze społeczeństwem w zakresie informowania policji o zaobserwowanych zagrożeniach w danym rejonie, w tym także tych związanych z przestępstwami przeciwko środowisku (identyfikowane jako dzikie wysypiska śmieci – patrz poniższa mapa).

Uprzednio informacje o rezultatach działania policji w latach 1999-2020 w zakresie przestępstw środowiskowych zebrałem w:

Obecnie moją intencją jest, zanim upublicznione zostaną statystyki za 2021 r., przedstawienie sposobu na prognozowanie ilości stwierdzonych przestępstw w gospodarce odpadami (na podstawie art. 183 kk) z wykorzystaniem klasycznej metody opartej na technice autoregresji i średniej ruchomej dla szeregu czasowego.

ARIMA – autoregressive integrated moving average

Mowa o metodzie określonej w języku angielskim jako autoregressive integrated moving average, w skrócie ARIMA, a tłumacząc na język polski o metodzie autoregresyjnej średniej ruchomej z różnicowaniem. Zatem już sama nazwa wskazuje, że mamy tu do czynienia z łączeniem kilku elementów: procesu autoregresji, średniej ruchomej oraz różnicowania, które prowadzi do przekształcenie szeregu niestacjonarnego w stacjonarny, co jest niezbędne dla prawidłowości funkcjonowania opracowywanego modelu. Ale po kolei:

Autoregresja

autoregresja to proces, w którym każda wartość jest liniową kombinacją poprzednich wartości, czyli korzysta niejako z pamięci procesu. Autoregresja może być różnego rzędu, a on określa ile wcześniejszych wartości ma wpływ na wartość bieżącą. Np. gdyby ujmować wpływ jedynie poprzedniej wartości proces autoregresyjny można by zapisać jako:

gdzie: y_t– wartość w chwili t, y_t-1 – wartość w chwili t-1, e_t – składnik losowy w chwili t

Wówczas oznacza się go jako AR(1). Jeżeli rząd wynosiłby 2 to proces AR(2) przedstawiałby się następująco:

Parametry j określają wpływ kolejnych, coraz bardziej oddalonych w szeregu wartości, na wartość bieżącą. W modelu ARIMA rząd procesu autoregresyjnego określany jest często literą p.

Średnia ruchoma

średnia ruchoma – jest doskonałym narzędziem śledzenia trendu z jednoczesnym odfiltrowaniem zbędnego szumu. W pewien sposób jest podobna do autoregresji, gdyż każda wartość elementu średniej ruchomej w modelu stanowi określoną q-elementową kombinację poprzednich wartości składnika losowego. Proces ten możemy zapisać jako:

gdzie: y_t– wartość w chwili t, e_t-1 – składnik losowy w chwili t-1, e_t – składnik losowy w chwili t i oznaczyć jako MA(q). Parametry b określają wpływ kolejnych, coraz bardziej oddalonych zaburzeń w szeregu na wartość bieżącą szeregu. Czyli możemy stwierdzić, że proces średniej ruchomej MA(q) jest średnią ważoną z ostatnich zaburzeń składników losowych modelu.

Proces integracji

różnicowanie – proces polegający na przekształceniu szeregu niestacjonarnego w stacjonarny, czyli taki, którego własności nie zmieniają się w czasie. Często możemy kojarzyć taki szereg jako taki, z którego wyeliminowano trend i sezonowość. W rezultacie szereg stacjonarny prezentuje się za pomocą stałej średniej i wariancji w czasie. Aby doprowadzić szereg niestacjonarny do stacjonarnego wymagane może być przeprowadzenie kilkukrotnego różnicowania, d.

Stąd łącząc wszystkie trzy składowe otrzymujemy model ARIMA o parametrach p, q i d, oznaczany jako ARIMA(p, d, q), który możemy przedstawić w postaci następującego równania:

gdzie y’ oznacza szereg czasowy różnicowany jednokrotnie. Jak wspomniałem różnicowanie prowadzimy do uzyskania stacjonarności szeregu. Dobór wartości parametrów p, q i d wyjaśnię na etapie realizowanego przykładu.

Model prognozowania ilości stwierdzonych przestępstw w gospodarce odpadami z art. 183 kk

Dane do modelowania pobrane zostały z policyjnego portalu Statystyka. Dane dostępne są w postaci arkusza Excel, i zawierają informacje o liczbie postępowań wszczętych, liczbie przestępstw stwierdzonych, liczbie przestępstw wykrytych oraz wskaźniku wykrycia.

Szereg czasowy liczby stwierdzonych przestępstw z artykułu 183 kodeksu karnego przedstawia się następująco:

A jeżeli przyjrzymy się typowi danych w poszczególnych kolumnach to mamy do czynienia w obu przypadkach z liczbami typu całkowite (int64). Przy okazji zerknijmy także na graficzną reprezentacje danych. Wyraźnie rzuca się w oczy ogromny wzrost liczby stwierdzonych przestępstw w latach 2018-2020. Jest to także przyczyna konieczności zaostrzenia kar za niezgodne z prawem działania w gospodarce odpadami poprzez wprowadzane obecnie przepisy ustawy o zmianie niektórych ustaw w celu przeciwdziałania przestępczości środowiskowej.

W modelowaniu szeregów czasowych istotny jest element czasu. Stąd tez kolumnie „Rok” zostanie nadany właściwy typ danych i zostanie ona zamieniona na kolumnę indeksową. Moduł pandas udostępnia metodę to_datetime (pandas.to_datetime – documentation) która sprytnie to realizuje. Przy okazji określiłem częstotliwość danych jako roczne podane na ostatni dzień grudnia.

Mamy zatem przygotowany szereg czasowy i czas zająć się określeniem parametrów p, d, q w modelu ARIMA.

Jak określić parametry w modelu ARIMA?

Zaczynijmy od parametru d. Aby go określić musimy ocenić stacjonarność szeregu. Do tego celu skorzystam z testu statystycznego ADF (ang. Augmented Dickey Fuller Test).

Procedura testu ADF, zwanego także testem pierwiastka jednorodnego, testuje hipotezę zerową (H₀: sig = 0) zakładającą, że proces ma pierwiastek stopnia pierwszego, natomiast hipoteza alternatywna (H₁: sig < 0) zakłada stacjonarność analizowanego procesu. W module statsmodels.tsa.stattools mamy metodę adfuller do realizacji testu ADF.

weryfikacja stacjonarności szeregu czasowego

Szereg stał się stacjonarny dopiero po sześciokrotnym różnicowaniu. Stąd też d = 6. Poza konieczności kilkukrotnego powtarzania operacji różnicowania czynność wyznaczenia parametru d była łatwa. Jeżeli chodzi o parametry p i q sprawa jest nieco trudniejsza, i wymaga abyśmy pomogli sobie funkcjami autokorelacji (ACF) oraz cząstkowej autokorelacji (PACF). Obie te funkcje znajdziemy w module statsmodels.tsa.stattools i ich znaczenie w identyfikacji parametrów p i q dobrze przedstawione zostało w artykule.

Aby ustalić wartość parametru p w modelu autoregresji AR(p) korzystamy z funkcji autokorelacji cząstkowej, gdyż można wykazać, że dla procesu AR(p) funkcja cząstkowej autokorelacji znika dla opóźnień większych niż rząd p. Zatem jeżeli cząstkowa autokorelacji PACF(h) znajduje się pomiędzy przedziałami ufności 1.96/sqrt(n) dla h>p to spodziewany się, że dane są realizacją procesu autoregresji rzędu p. Stąd p = 1.

Podobnie jeżeli autokorelacji ACF(h) znajduje się pomiędzy przedziałami ufności 1.96/sqrt(n) dla h>p to spodziewany się, że dane są realizacją procesu średniej ruchomej q. Stąd q = 0.

Zatem parametry końcowe modelu są następujące: p = 1, d = 6, q = 0 tworząc model ARIMA(1,6,0).

Dopasowanie modelu ARIMA

Do dopasowania modelu do danych szeregu skorzystam z klasy ARIMA z modułu statsmodels.tsa.arima_model przypisując wyznaczone hiperparametry p, d i q. Wywołanie metody fit tego modułu zwraca nam model, który prezentuje się następująco:

W wyniku uzyskujemy stałą i współczynniki do równania opisującego model. W związku z tym, ze określiliśmy parametr q na 0 nie mamy składnika średniej ruchomej, natomiast jeden współczynnik (p = 1) identyfikuje proces autoregresji. Ponadto model dostarcza nam kilku statystycznych miar do jego oceny. Wskaźnik AIC (ang. Akaike Information Critera) uważa się za dobry miernik matematycznego dopasowania i powinien być możliwie najmniejszy. Klasa ARIMA umożliwia podanie parametrów p i q w postaci listy wartości celem znalezienia optymalnego rozwiązania. Mamy też możliwości optymalnego wyboru modelu. W pakiecie pmdarima opracowano metodę automatycznego doboru parametrów. W oparciu o optymalizator Broyden-Fletcher-Goldfarb-Shanno metoda potwierdziła, że optymalny model to ARIMA(1,6,0).

Jeżeli spojrzy się na oryginalny wykres pełnego szeregu czasowego to łatwo zauważyć wykładniczy wzrost ilości stwierdzonych przestępstw w okresie ostatnich lat, czyli 2019-2020. To jest pewnego rodzaju problem, model ARIMA ma trudne zadanie do przewidzenia, gdyż jest to coś nietypowego w niewielkim zbiorze danych. Zresztą widać to po uzyskanych wartościach prognozowanych.

Stąd zrezygnowałem nawet z obliczania statystycznych metryk porównujących dane rzeczywiste z wartościami otrzymanymi na podstawie modelu prognostycznego.

Chciałoby się rzec, że gdyby nie lata 2019-2020 to tworzony model prezentowałby się zupełnie inaczej. A tak, model oparty na liniowej zależności, taki jak ARIMA, nie spełnia swojego zadania w tym miejscu. Zasadniczo skuteczniejsze i poprawne byłoby wykorzystanie modelu nieliniowego, np poprzez zastosowanie modelu opartego na sieciach neuronowych LSTM. Będzie to ciekawy materiał na następny wpis.

Tymczasem udając się na wakacyjny odpoczynek 😊 życzę sobie i wszystkim, którzy wytrwali do końca tego wpisu udanej aury pogodowej. Być może w między czasie statystyki policyjne za rok 2021 zostaną zaktualizowane i będziemy prognozować wartości na lata kolejne.