Ekonometryczne modelowanie popytu na energię elektryczną ..., Streszczenia z Transport

Pobierz Ekonometryczne modelowanie popytu na energię elektryczną ... i więcej Streszczenia w PDF z Transport tylko na Docsity! Uniwersytet Warszawski Wydział Nauk Ekonomicznych Rafał Woźniak Ekonometryczne modelowanie popytu na energię elektryczną za pomocą danych wysokiej częstotliwości Praca doktorska Praca wykonana pod kierunkiem dr. hab. Ryszarda Kokoszczyńskiego, prof. UW z Zakładu Finansów Ilościowych WNE UW Warszawa, czerwiec 2015 2 Podziękowania Dziękuję prof. Ryszardowi Kokoszczyńskiemu za pomoc w nadaniu ostatecznej struktury pracy oraz za cenne uwagi. Dłużny jestem podziękowania dr Natalii Nehrebeckiej za cier- pliwe odpowiedzi na moje niekończące się pytania. Marcinowi Sienickiemu dziękuję za owocną dyskusję na temat zmiennych użytych w modelowaniu wskaźników koniunktury. Pani Annie Włostowskiej i Krystynie Karpiuk dziękuję za zrozumienie. Słowa podzięko- wania należą się Wojciechowi Broniarkowi za bycie pierwszym krytycznym czytelnikiem opisów. Szczególnie dziękuję Marcie za motywację. Wstęp Niewątpliwie energia elektryczna we współczesnym świecie odgrywa ogromną rolę. Specy- fika elektryczności, a w zasadzie brak możliwości jej ekonomicznego magazynowania, spra- wia, że analizowanie popytu na energię elektryczną jest niezmiernie ważne dla wszystkich podmiotów z rynku energii. Na rynku energii elektrycznej w Polsce, od momentu jego liberalizacji 1 lipca 2007, wyróżnia się wytwórców energii, pośredników, firmy zajmujące się przesyłem energii oraz odbiorców. Dysponowanie idealnymi prognozami zapotrze- bowania na energię elektryczną mogłoby pozwolić każdemu uczestnikowi tego rynku na ekonomiczne gospodarowanie, ponieważ każde niedostateczne zakontraktowanie energii lub przeszacowanie zapotrzebowania generuje koszty. W centrum zainteresowania są wy- sokoczęstotliwościowe charakterystyki rynku energii, a dane 15-minutowe wydają się być idealne do tego typu analiz. W niniejszej pracy sprawdzone zostały dwie główne hipotezy badawcze. Pierwsza z nich dotyczyła istnienia nieliniowego wpływu temperatury powietrza na zapotrzebowanie na energię elektryczną w obserwowanych danych 15-minutowych. Drugą zasadniczą hipotezą badawczą jest istotność pozytywnego wpływu cyklu koniunkturalnego na analizowane za- potrzebowanie na moce przesyłowe krajowych sieci elektroenergetycznych czyli sensow- ność wykorzystania prognoz zagregowanego stanu gospodarki do sporządzania prognoz średniookresowych. Wybór tych dwóch hipotez uargumentowany jest chęcią formalnego sprawdzenia postaci i pokazania nieliniowej zależności między temperaturą a zapotrze- bowaniem na energię, natomiast druga, związana z cyklem koniunkturalnym, wydaje się niezwykle ciekawa dla autora. Poza wymienionymi hipotezami głównymi wartymi rozważenia wydają się być jeszcze do- datkowe kwestie. Pierwszą hipotezą poboczną jest istnienie w obserwowanych danych kwadransowych opisujących zapotrzebowanie na moce przesyłowe potrójnej struktury se- zonowej związanej z wahaniami rocznymi, tygodniowymi oraz fluktuacjami charaktery- stycznymi dla poszczególnych godzin. Drugim wartym weryfikacji stwierdzeniem jest to o istnieniu w rozważanych danych tendencji wzrostowej - trendu. Trzecim dodatkowym zagadnieniem rozważanym w pracy jest weryfikacja przypuszczenia o tym, że metoda prognozowania oparta o modele regresji liniowej dla poszczególnych kwadransów powinna okazać się najlepsza dla polskich danych. 5 6 SPIS TREŚCI Hipotezy główne i dodatkowe wydają się być ekonomicznie intuicyjne. Obserwowany wzrost gospodarczy Polski w okresie 2002-20141 przekładać się musi na wzrost zapotrzebo- wania na energię elektryczną. Oczekiwanie skomplikowanej struktury wahań okresowych zdaje się być bezdyskusyjne dla każdego analityka i uczestnika rynku energii. Cykle dni i nocy, tygodniowe struktury dni pracy i weekendów oraz serie pór roku, czyli zmiennej dłu- gości dnia, uzasadniają postawienie hipotezy o potrójnej strukturze wahań sezonowych. Hipoteza o przekładaniu się koniunktury na potrzeby elektroenergetyczne gospodaruja- cych podmiotów wydaje się mieć uzasadnienie w teorii mikroekonomii. Przedsiębiorca lub firma planuje produkcję w zależności od bieżącego i spodziewanego popytu. Słabszy zagregowany popyt związany jest z okresami gorszej koniunktury, więc rozsądnym wydaje się oczekiwanie przekładania się sytuacji gospodarczej na energię elektryczną. W niektó- rych przypadkach wydawać się może, że zależność ta ma przeciwny zwrot, gdyż proces produkcyjny wyprzedza moment sprzedaży. Dokładnie z tego samego powodu oczekiwać można logiczności wykorzystania prognoz sytuacji ekonomicznej do wyznaczania przy- szłych potrzeb energetycznych podmiotów. Zdecydowana większość badań wpływu cyklu koniunkturalnego na popyt na energię elektryczną wykorzystuje dane roczne i kwartalne. Wartość dodana niniejszej pracy polega także na analizowaniu wspomnianej relacji dla danych o wyższej, miesięcznej częstotliwości. W niniejszej pracy podstawą analizy są dane 15-minutowego zapotrzebowania na moce przesyłowe krajowych sieci elektroenergetycznych (KSE). Informacje te bez wątpienia można określić mianem wysokoczęstotliwościowych. Dane są dostępne dla całej sieci, więc poziomem zagregowania jest cały kraj. Zamiennie będą używane pojęcia zapotrzebowania na moce przesyłowe KSE, popytu na energię i popytu na moce przesyłowe. Definicyjnie nie są to tożsame pojęcia, jednak ze względu na poziom zagregowania danych (cały kraj) pojęcia te mogą być postrzegane jako tożsame. Pomijana jest tutaj kwestia eksportu i importu energii, ponieważ w wykonanych analizach nie dokonywane było rozróżnienie pochodzenia energii w sieci. Praca jest skonstruowana tak, że kolejno opisywane są poszczególne składowe szeregu 15-minutowego zapotrzebowania na moce przesyłowe KSE i weryfikowane kolejne hipo- tezy. Pierwszy rozdział zawiera krótki opis analizowanej zmiennej oraz jej podstawowe charakterystyki. Głowna część tego rozdziału poświęcona jest sprawdzeniu istnienia ten- dencji wzrostowej oraz analizom stacjonarności szeregu. Dla tych drugich zaproponowana została specjalna wersja testu stacjonarności sezonowych danych oraz zaprezentowane zo- stały jego stablicowane wartości krytyczne. Wykorzystane zostało ponadto podstawowe narzędzie ekstrakcji komponentu trendu z szeregu oraz przedstawione zostały jego para- metry, jakie powinny zostać zastosowane w przypadku danych 15-minutowych. W rozdziale drugim zaprezentowane zostały analizy struktury sezonowej szeregu zapo- 1Silna tendencja wzrostowa w polskim produkcie krajowym brutto, związana z transformacją gospo- darki, rozpoczęła się w połowie lat 90’, jednak ze względu na zakres czasowy analizowanych danych przytoczony jest krótszy okres. SPIS TREŚCI 7 trzebowania na moce przesyłowe krajowych sieci elektroenergetycznych. Nie tylko spraw- dzone zostało istnienie wielorakiej struktury wahań okresowych, ale także przeanalizowane zostały wersje tych wahań. Wychodząc od podstawowego modelu wygładzania wykładni- czego zaprezentowane zostały jego kolejne rozszerzone wersje wyestymowane dla polskich danych. Zaprezentowane zostało dodatkowe podejście, które, jak się wydaje, nie zostało jeszcze wykorzystane w literaturze. Rozdział trzeci poświęcony jest analizom cyklu koniunkturalnego i podstawowym sposo- bom sprawdzenia jego przekładania się na popyt na energię elektryczną. Ze względu na to, że w literaturze polskiej nie istnieje wiele opracowań zależności między zagregowanym sta- nem gospodarki a popytem na energię elektryczną w sieci zaprezentowane zostały powody, dla których istotne jest rozważanie tej kwestii oraz uwzględnienie jej w procesie mode- lowania i przede wszystkim prognozowania. Zaprezentowany został sposób oszacowania indeksu przybliżającego cykl koniunkturalny oraz jego wersja pozwalająca prognozować te wartości. Rozdział zawiera także podrozdział poświęcony propozycji szacowania sze- regu miesięcznego produktu krajowego brutto, który mógłby najdokładniej odzwierciedlać sytuację w gospodarce. Wykorzystanie zmiennych egzogenicznych w modelowaniu zapotrzebowania na energię elektryczną jest tematyką rozdziału czwartego. Zaprezentowane i wykorzystane zostały stochastyczne i deterministyczne zmienne objaśniające, czyli zmienne pogodowe i dni wolne od pracy. Ponownie, ze względu na nieliczność badań tego typu dla Polski, wyko- nana została szeroka analiza dla częstotliwości miesięcznych, tygodniowych i dziennych. Sprawdzone zostało istnienie nieliniowej zależności między temperaturą a energią za po- mocą modeli parametrycznych i nieparametrycznych. Rozdział piąty zawiera eksperymenty, w których porównane zostały sposoby prognozowa- nia zapotrzebowania na energię. Rozważone zostały dwa horyzonty prognoz: 3-dniowy i 3-miesięczny. Wybór okresu 3-dniowego jest uargumentowany ogromnym zainteresowa- niem teoretyków i praktyków prognozami w tym horyzoncie, natomiast okres 3-miesięczny umożliwia weryfikację hipotezy o użyteczności informacji o cyklu koniunkturalnym do prognozowania średniookresowego. Wykorzystane zostały modele uznane w literaturze za najdokładniejsze dla najczęściej opisywanego w literaturze rynku energii elektrycznej w Stanach Zjednoczonych. Kończące pracę zakończenie zawiera rekapitulację najważniejszych wniosków i wyniki we- ryfikacji hipotez. Przedstawione zostały także kwestie wymagające dalszych prac i analiz w tej tematyce oraz potencjalne rozszerzenia zaprezentowanych metod. Praca ma charakter akademicki z pewnymi elementami praktycznymi. Praktycy mogą mieć wypracowane własne metody i sposoby, które pozwalają uzyskać bardzo dobre rezul- taty. Osoby zajmujące się zawodowo prognozowaniem popytu na energię elektryczną mogą mieć dostęp do komputerów z większymi możliwościami obliczeniowymi, dzięki czemu 10 ROZDZIAŁ 1. ANALIZA JEDNOZMIENNOWA dostępnej za 1 stycznia 2002 r. do 10 października 2013 roku. Dane są danymi wysoko- częstotliwościowymi, ponieważ są próbkowane co kwadrans. W związku z tym, w próbie znajdują się 412892 obserwacje. Analiza graficzna szeregu, przedstawionego na rysunku 1.1, wskazuje na występowanie tendencji wzrostowej w analizowanym 12-letnim okresie. Ze względu na przejrzystość na rysunku 1.1 zaprezentowane zostały tylko dane dla trzecich śród każdego miesiąca. Uzasadnione może być to ponadto tym, że zgodnie z nieformalną regułą postępowania za- pewnić ma to lepszą porównywalność poprzez usunięcie części efektów sezonowych zabu- rzającyh wnioskowanie. Zauważalne jest także zróżnicowanie wartości dla poszczególnych godzin oraz cykle większego i mniejszego zapotrzebowania zimą i latem. Wydaje się także, że istotny spadek miał miejsce w 2009 roku. Powyższe obserwacje zostały zweryfikowane w niniejszym, drugim i trzecim rozdziale. Zapotrzebowanie na moce przesyłowe w krajowym systemie elektroenergetycznym jest zmienną stanu, w przeciwieństwie do wielu kategorii makroekonomicznych, które są prze- pływami. Oznacza to, że szereg ten musi być analizowany podobnie jak temperatura, stan wody w rzece czy kurs walutowy. Zmiana częstotliwości na niższą nie odbywa się, jak dla przepływów, poprzez sumowanie, ale poprzez wybranie jednej z obserwacji z analizowa- nego okresu, na przykład średniego poziomu lub wartości na koniec okresu. W niniejszej pracy zajmuję się modelowaniem zapotrzebowania na moce przesyłowe kra- jowego systemu energetycznego. Niezwykle ważne dla wielu uczestników rynku energii jest także modelowanie i prognozowanie cen energii elektrycznej. Obszerną analizę metod prognozowania cen energii przedstawił Weron2. 1.2 Analiza spektralna Analiza częstotliwościowa pozwala zbadać występowanie wahań o poszczególnych często- ściach w szeregu czasowym. Dane o zapotrzebowaniu na moce krajowej sieci elektroener- getycznej potraktowane zostały w poniższych analizach jako dane dzienne, próbkowane 96 razy na każdy dzień, czyli co kwadrans. Alternatywnie uznać można, że okres wynosi tydzień, a próbkowanie równa się 672 razom w tygodniu. Wybór okresu nie wpływa na wnioski, a jedynie na sposób interpretacji wyników. W typowych analizach zmiennych miesięcznych lub kwartalnych za okres przyjmuje się odpowiednio 12 lub 4 obserwacje, natomiast w przypadku danych o złożonej strukturze sezonowej wybór ten nie jest jedno- znaczny. W widmie szeregu, przeskalowanego do skali logarytmicznej, widoczne są składowe od- powiadające wahaniom długookresowym. Potwierdza to tezę o niestacjonarności szeregu. 2Weron R., 2014, Electricity Price Forecasting: A Review of the State-of-the-Art with a Look into the Future, International Journal of Forecasting, 30, s. 1030-1081. 1.2. ANALIZA SPEKTRALNA 11 Rysunek 1.2: Widmo oryginalnego szeregu - skala logarytmiczna Źródło: Opracowanie własne. Argumenty widma rozważanego szeregu przebiegają przedział [0, 48], gdzie 48 jest często- tliwością Nyquista. Transformacji Fouriera poddawane są szeregi stacjonarne, w związku z czym z danych usunięty został deterministyczny trend liniowy. Uzasadnione jest to tym, że takie działanie nie zmienia charakterystyki częstotliwościowej szeregu poza nieskończo- nymi częstotliwościami. Pik w widmie oszacowanym dla częstości 1 odpowiadającej wahaniom powtarzającym się 1 raz na 96 kwadransów świadczy o występowaniu sezonowości przypisywanej poszczegól- nym kwadransom, czyli sezonowości dziennej. Piki powtarzające się dla kolejnych całko- witoliczbowych częstości odpowiadają harmonicznym sezonowości dziennej i wynikają z własności transformaty Fouriera. Zawężając analizę gęstości spektralnej procesu do przedziału [0, 1], wnioskowanie zostaje ograniczona do wahań rzadszych niż 1 raz na dzień, co prezentuje rysunek 1.4. Zauważyć można lokalne maksimum mocy dla argumentu 0,1428, co jest przybliżeniem 1 7 , czyli wahaniom powielającym się co 7 dni i związane z sezonowością tygodniową. Piki dla częstości 2 7 , 3 7 , . . . , 6 7 są harmonicznymi tej sezonowości. Rysunek 1.5 przedstawia widmo ograniczone do jeszcze węższego przedziału. Lokalne maksimum funkcji mocy dla częstości 0,00279, oraz jej harmonicznych, związane jest z wahaniami powielającymi się co około 365,25 dnia. Poprawka o 0,25 wynika z kwestii roku przestępnego pojawiającego się w próbce raz na 4 lata. Podejście to jest stosowane przez autorów zajmujących się tematyką sezonowości w danych dziennych, na przykład Dordonnat i inni3. Warto zauważyć także wysokie wartości widma dla niskich częstości. Wydaje się, że - 3Dordonnat V., Koopman S.J., Ooms M., Dessertaine A., 2008, An Hourly Periodic State Space Model for Modeling French National Electricity Load, International Journal of Forecasting, 24, 588-587. 12 ROZDZIAŁ 1. ANALIZA JEDNOZMIENNOWA Rysunek 1.3: Widmo szeregu z odjętym trendem liniowym - skala logarytmiczna Źródło: Opracowanie własne. oprócz trendu - w danych występować może składnik cykliczny, odpowiadający waha- niom średniookresowym. Na rysunku 1.4 świadczyć o tym może lokalne maksimum mocy dla częstości 0,001163, czyli dla wahań co 2,35 roku. Zgodnie z najczęściej stosowaną definicją cyklu koniukturalnego Burnsa i Mitchella4 za wahania odpowiadające cyklowi koniunkturalnemu uznać należy te, które powtarzają się w zakresie od półtora roku do 8 lat. Uzyskana wartość wpada w definicyjny przedział. Obecność zaburzeń w obserwowanym szeregu może wpływać na uzyskiwane widmo mocy. Powoduje to problem tzw. wycieku mocy lub ucieczki częstotliwości (spectral leakage), który skutkuje przeniesiem części mocy z właściwej częstotliwości np.: ω = ω0 do innej ω = ωa. Rozwiązaniem tego problemu jest zastosowanie okienkowania sygnału. W świetle twierdzenia o splocie przemnożenie oryginalnego sygnału przez okienko jest równoznaczne w domenie częstotliwości ze splotem jego funkcji mocy z okienkiem. Problem wycieku mocy został opisany między innymi przez Harrisa5, który też porównał różne okienka. Sygnał został zokienkowany trzema różnymi, często wykorzystywanymi, oknami Blackmana- Harrisa, Hamminga oraz Tuckeya. Uzyskana w ten sposób gęstość widmowa sygnału jest bardzo podobna w każdym z tych przypadków do oryginalnego sygnału. Rysunek 1.6 prezentuje część widm mocy uzyskanych z zastosowania trzech wspomnianych okien. W przypadku okien Blackmana-Harrisa i Hamminga obserwowany jest pik dla częstości 0,001163, odpowiadającej wahaniom powtarzającym się co 2, 35, co daje podstawę do dal- szej analizy obecności wahań średniookresowych w szeregu zapotrzebowania mocy krajo- 4Burns A.F., Mitchell W.C., 1938, Statistical Indicators of Cyclical Revival, National Bureau of Eco- nomic Research, Nowy Jork. 5Harris , 1976, Windows, Harmonic Analysis, and the Discrete Fourier Transform, Naval Undersea Center, San Diego. 1.3. STACJONARNOŚĆ I SEZONOWOŚĆ 15 Tablica 1.1: Modele testujące stacjonarność i sezonowość danych (1) (2) (3) (4) zmienna ln (yt) ln (ybtt ) ∆yt ∆yt stacjonarność δ - - -53,51 -24,59 stała tak tak tak nie opóźnienia zmiennej zależnej 1;5 1;5 1/35 1/35 test braku sezonowości statystyka testowa 348,61 348,18 - - p-value 0,000 0,000 - - test braku autokorelacji statystyka testowa 8,425 9,571 1,018 3,197 p-value 0,004 0,002 0,313 0,074 Źródło: Opracowanie własne. Pozostałe kolumny tablicy 1.1 poświęcone są wynikom badania rzędu stacjonarności. Naj- częściej stosowanym testem do testowania stopnia integracji jest test Dickey’a-Fullera7, w skrócie DF. Testowana jest hipoteza H0 : δ = 0, przeciwko alternatywie H1 : δ < 0 w równaniu ∆yt = δyt−1 + εt. (1.1) Test nazywany jest też testem pierwiastka jednostkowego, powieważ równanie (1.1) może zostać zapisane w postaci yt = ρyt−1 + εt, (1.2) gdzie ρ = (1 + δ), która jest równoważna 1.1. Odrzucenie hipotezy zerowej H0 na ko- rzyść alternatywy oznacza, że szereg yt jest zintegrowany w stopniu zero. Ze względu na obciążenie estymatora parametru δ sposób weryfikacji założonej hipotezy jest odmienny od standardowo przyjętego w modelach regresji liniowej. Rozkład statystyki testowej t = δ̂/se(δ̂) nie jest rozkładem t-Studenta, a jego tablice mogą być generowane za po- mocą metody Monte Carlo. Tablice rozkładu statystyki testowej zosatły przedstawione przykładowo przez Charemzę i Deadmana8. W przypadku sezonowego szeregu czasowego w równaniu (1.1) pojawiają się dodatkowe regresory odpowiadające sezonowości. Dla danych 15-minutowych z wahaniami sezonowymi równanie (1.1) przyjmie postać: ∆yt = δyt−1 + 96∑ i=1 QK,i,t + 6∑ i=1 QD,i,t + 11∑ i=1 QM,i,t + εt, (1.3) 7Dickey D.A., Fuller W.A., 1979, Distribution of Estimators for Autoregressive Time Series With a Unit Root, Journal of American Statistical Association, Vol. 74, No. 366, s. 427-431 8Charemza W., Deadman D. 1997, Nowa ekonometria, Polskie Wydawnictwo Ekonomiczne, s. 292 16 ROZDZIAŁ 1. ANALIZA JEDNOZMIENNOWA gdzie QK,i,t, QD,i,t, QM,i,t są zmiennymi zero-jedynkowymi odpowiadającymi odpowiednio poszczególnym kwadransom, dniom i miesiącom. Na podstawie modeli (3)-(4) wyciągnąć należy wniosek o zerowym stopniu zintegrowania szeregu. Konkluzja ta jest sprzeczna z analizą graficzną szeregu i przypuszczeniami. Prawdopodobną przyczyną takiego rezul- tatu może być występowanie autokorelacji wysokiego rzędu, występowanie także stocha- stycznej wersji sezonowości, stochastycznego trendu lub niska moc testu DF, a rozwiąza- niem umożliwiającym doprecyzowanie wniosków może być bardziej skomplikowany test pierwiastka jednostkowego. W związku z wymienionymi słabościami nie można poprzestać na tym podejściu należy wykorzystać inne narzędzia. Test zaproponowany przez Hylleberga i innych9, zwany da- lej HEGY, jest testem pozwalającym zidentyfikować rodzaj sezonowości i jednocześnie zweryfikować hipotezę o istnieniu pierwiastka jednostkowego. Test HEGY dla danych 15-minutowych opiera się na modelu regresji postaci: ∆96yt = 96∑ i=1 aiQi,t + 96∑ i=1 biYi,t−1 + k∑ i=1 ci∆96yt−i + εt, (1.4) gdzie k jest liczbą opóźnień zmiennej objaśnianej, zmienne Qi są sezonowymi zmiennymi zero-jedynkowymi, natomiast Yi,t są zmiennymi zdefiniowanymi w następujący sposób: Y1,t = 96∑ j=1 cos (j · 0)yt−j, (1.5) Y2i,t = 96∑ j=1 cos (j2π · (2i)/96)yt−j, (1.6) Y2i+1,t = 96∑ j=1 sin (j2π · (2i+ 1)/96)yt−j, (1.7) dla i = 1, 2, . . . , 47 oraz Y96,t = 96∑ j=1 cos (j2π · 96/96)yt−j. (1.8) Zgodnie z interpretacją Charemzy i Deadmana10, gdy wszystkie współczynniki bi są równe zeru, a współczynniki ai jednakowe, to w obserwowanym szeregu yt występuje tylko sto- chastyczna sezonowość. W przypadku różnych wartości ai i przynajmniej jednego bi różnego od zera obserwowana jest zarówno stochastyczna, jak i deterministyczna wer- sja wahań sezonowych. Poszczególne zmienne Yi odpowiadają wahaniom z określonymi częstotliwościami. Tak na przykład Y1 odpowiada integracji szeregu, a obydwie Y2 oraz Y3 wahaniom sezonowym powtarzającym się co 96 okresów 15-minutowych. 9Hylleberg S., Engle R.F., Granger C.W.J., Yoo B.S., 1990, Seasonal Integration and Cointegration, Journal of Econometrics, Vol. 44(1-2), s. 215-238 10Charemza W., Deadman D. 1997, Nowa ekonometria... 1.4. ANALIZA KOMPONENTU TRENDU 17 Tablica 1.2: Statystyki testowe i wartości krytyczne testu HEGY Wartość Hipoteza ∆96yt krytyczna 5% H1 0 : b1 = 0 28,625 -2,82 H2 0 : b2 = b3 = 0 1654,452 3,37 H3 0 : ∀i ai = 0 26278,304 3,43 Źródło: Opracowanie własne. Z powodu braku programu umożliwaiającego przeprowadzenie testu HEGY dla danych 15-minutowych konieczne było samodzielne oprogramowanie tego testu11. Tablica 1.2 prezentuje wyniki testu HEGY dla badania stacjonarności analizowanych danych 15- minutowych. Dodatnia wartość statystyki testowej w testowaniu H1 0 pozwala stwierdzić występowania pierwiastka jednostkowego, co potwierdza dodatkową hipotezę badawczą o występowaniu rosnącego trendu. Potwierdzić należy występowanie stochastycznej se- zonowości ze względu na wartość statystyki testowej 1654,452 wykraczającej poza 95% przedział ufności. Co więcej, uznać należy także fakt istotności deterministycznej sezono- wości w rozważanym szeregu. Test HEGY został przeprowadzony raz jeszcze, z tą różnicą, że zmienną objaśnianą w rów- naniu (1.4) były różnice ze skorygowanym odjemnikiem ∆̃96yt. Różnice ze skorygowanym odjemnikiem różnią się od klasycznych różnic rzędu 96 tym, że dla okresów związanych ze zmianą czasu odjemnikiem jest 95 lub 97 opóźnienie zmiennej. Korekta ta służyć ma sprawdzeniu czy zmiana czasu związana z wprowadzaniem czasu letniego nie zaburza wnioskowania. Wnioski z testu HEGY dla właściwego różnicowania zmiennej pozostają jednak niezmienne. 1.4 Analiza komponentu trendu Filtr Hodricka-Prescotta12 został zaproponowany około 18 lat temu i obecnie jest jednym ze standardowych narzędzi ekonometrycznych. Narzędzie to służy do wyodrębniania z szeregu czasowego długookresowego trendu. Różnica między oryginalnym i wynikowym szeregiem stanowi wahania cykliczne i nieregularne. Filtr Hodricka-Prescotta, dalej HP, jest liniowym, dwustronnym filtrem za pomocą któ- rego można wyznaczyć wygładzony komponent µ szeregu y. Szkielet filtru stanowi przed- stawienie szeregu czasowego yt w postaci sumy komponentu wzrostowego µt oraz kompo- 11Stworzony w tym celu został skrypt ”HEGY96.ado”, który jest instalowalną procedurą dla pakietu Stata. 12Hodrick R.J., Prescott E.C., 1997, Postwar U.S. Business Cycles: An Empirical Investigation, Jour- nal of Money, Credit and Banking, Vol. 29, No. 1, s. 1-16. 20 ROZDZIAŁ 1. ANALIZA JEDNOZMIENNOWA Rysunek 1.7: Szereg zapotrzebowania na moce przesyłowe i jego trend Źródło: Opracowanie własne. z powyższym należy zastosować regułę Ravna i Uhliga, w której λ = 1600 · ( l 4 )4 , gdzie l jest liczbą okresów w roku. Identyczny wynik otrzymuje się stosując wzór przed- stawiony przez Harvey’a i Trimbura20, gdzie λ = [ 2 sin (ωc 2 )]4 , gdzie ωc jest wybraną częstością odcięcia. Estymacja komponentu trendu µ została wykonana w programie MatLab na podsta- wie własnej implementacji modelu. W zastosowanym filtrze wartość l była równa l = 4 · 24 · 365 · 10, co odpowiada ωc = 2π 4·24·365·10 . Funkcja wiarygodności została zmaksy- malizowana za pomocą algorytmu wyżarzania symulowanego, a uzyskana wartość funkcji log-wiarygodności równa była 161310,60. Oszacowana wartość wariancji składnika loso- wego σ̂2 varepsilon =0,026798, przy odchyleniu standardowym tego oszacowania 4.0251·10−9. Statystyka istotności oszacowania jest argumentem potwierdzającym istnienie niestacjo- narnego komponentu w danych o zapotrzebowaniu na moce krajowej sieci energetycznej. Dodatkowo, rysunek 2 potwierdza założoną hipotezę o występowaniu pierwiastka jednost- kowego w analizowanym szeregu. 20Harvey A.C., Trimbur T.M., 2003, General Model-Based Filters for Extracting Cycles and Trends in Economic Time Series, The Review of Economics and Statistics, 85(2), str 245. 1.4. ANALIZA KOMPONENTU TRENDU 21 1.4.1 Zmiana czasu i brakujące obserwacje Zgodnie z Dyrektywą 2000/84/WE Parlamentu Europejskiego i Rady21 z dnia 19 stycz- nia 2001 r. w sprawie ustaleń dotyczących czasu letniego w krajach należących do Unii Europejskiej, w każdą ostatnią niedzielę marca następuje zmiana czasu uniwersalnego o godzinę, czyli przejście z czasu zimowego na letni. W związku z tym uznaje się, że po godzinie 1:00 następuje godzina 3:00, w efekcie doba tego dnia trwa 23 godziny. Uzasad- nione jest to lepszym dostosowaniem aktywności ludzi do godzin, w których jest najwięcej światła słonecznego. Zmiana czasu letniego na zimowy ma miejsce w ostatnią niedzielę października i polega przestawieniu zegarów o 60 minut do tyłu o godzinie 3:00. Aries i Newsham22 przytoczyli dane, które wskazują, że w krajach uprzemysłowionych ko- sumpcja elektryczności związana z oświetleniem odpowiada za 5-15% całkowitego zużycia energii elektrycznej23. Przywołany przez Aries i Newshama24 raport Instytutu Polityki Fi- skalnej w Indianie nie wskazał jednoznacznych wniosków dla tamtejszych danych, podczas gdy późniejsze badanie Kotchena i Grant25 zaprzeczyło istnieniu oszczędności wynikają- cych ze zmian czasu i wskazało na istotny wzrost zużycia energii tym wywoływany dla danych dla tego samego stanu. Badania przedstawione w Aries i Newsham26 oraz przez Hilla i innych27 wskazują na niejednoznaczność efektów wprowadzania czasu letniego. Ra- cjonalność wprowadzania czasu letniego nie może być analizowana tylko i wyłącznie na podstawie statystyk zużycia energii elektrycznej. Aries i Newsham odnotowali, że dodat- kowa godzina wieczorem przed zmierzchem w niektórych częściach Stanów Zjednoczonych przełożyła się na wzrost zużycia paliw poprzez zmianę zachowań ludzi. Wspomnieć tutaj także należy analizę Ferguson i innych28 oraz raport Królewskiego Towarzystwa Zapo- biegania Wypadkom przytoczony w Aries i Newsham29 pokazujące pozytywną zależność między czasem letnim a zmniejszeniem liczby ofiar wypadków samochodowych i potrą- ceń pieszych i rowerzystów. Dyskusja na temat zmiany czasu jest bardzo szeroka i nie wskazuje na jednoznaczne wnioski, więc weryfikacja racjonalności czasu letniego w Polsce pozostanie tematyką dalszych prac badawczych. Zagadnienie zmiany czasu nie było rozpatrywane w uprzednio przedstawionych modelach, 21Parlament Europejski i Rada Unii Europejskiej, 2001, Dyrektywa 2000/84/WE Parlamentu Europej- skiego i Rady z dnia 19 stycznia 2001 r. w sprawie ustaleń dotyczących czasu letniego, Dziennik Urzędowy Unii Europejskiej, 12 t.2. 22Aries M.B.C., Newsham G.R., 2008, Effect of Daylight Saving Time on Lighting Energy Use: A Literature Overview, Energy Policy, 36, s. 1858-1866. 23Aries M.B.C., Newsham G.R., 2008, Effect of Daylight Saving ..., str. 1858 24Tamże, str. 1862 25Kotchen M.J., Grant L.E., 2008, Does Daylight Saving Time Save Energy? Evidence from a Natural Experiment in Indiana, NBER Working Paper 14429. 26Aries M.B.C., Newsham G.R., 2008, Effect of Daylight Saving ..., str. 1858 27Hill S.I., Desobry F., Garnsey E.W., Chong Y.-F., 2010, The Impact on Energy Consumption of Daylight Saving Clock Changes, Energy Policy 38, s. 4955-4965. 28Ferguson S.A., Preusser D.F., Lund A.K., Zador P.L., Ulmer R.G., 1995, Daylight Saving Time and Motor Vehicle Crashes: The Reduction in Pedestrian and Vehicle Occupant Fatalities, American Journal of Public Health, Vol. 85, No. 1, s. 92-95. 29Aries M.B.C., Newsham G.R., 2008, Effect of Daylight Saving ..., str. 1863 22 ROZDZIAŁ 1. ANALIZA JEDNOZMIENNOWA ponieważ w modelach regresji liniowej każdemu kwadransowi przypisana była zmienna ze- rojedynkowa. Nie pojawiał się w nich problem następowania godziny 3:00 po 1:00. Inaczej jest w modelach szacowanych na szeregach czasowych. Jeśli model szacowany jest itera- cyjnie, tak jak na przykład modele przestrzeni stanów, wówczas może zachodzić potrzeba bardziej wyszukanego podejścia. Zilustrowaniu zagadnienia służyć może model autore- gresyjny rzędu 96 dla danych kwadransowych. W przypadku zmiany czasu prosty model AR(96) nie będzie poprawnie wyjaśniał zmienności modelowanej zmiennej. Wydaje się, że rozwiązaniem tego problemu mogłoby być potraktowanie czterech kwadransów pominię- tej godziny drugiej jako braki obserwacji. Podejście takie wymaga rozszerzenia modelu, które jest możliwe w elastycznych modelach estymowane za pomocą filtru Kalmana. Rozważmy podobnie jak w poprzedniej części pracy model określony równaniami (1.11)- (1.13), z tą różnicą, że zmienna yt zawiera braki obserwacji. Na skutek braku wartości zmiennej yt w niektórych okresach czasu niezbędne jest wprowadzenie zmian do modelu przestrzeni stanów i zmodyfikowanie filtru Kalmana. Niech y†t będzie zmienną zdefinio- waną dla wszystkich okresów czasu, w tym także dla pominiętych godzin drugich przy zmianie czasu, tak, że y†t = yt dla punktów czasu, gdy yt jest obserwowana i nieobser- wowaną w pozostałych przypadkach. W porównaniu do modelu (1.11)-(1.13) pojawia się dodatkowe równanie y†t = µt + εt, dla t = 1, . . . , T̃ . Dodatkowo, ze względu na brak możliwości obserwowania wartości zmiennej yt niezbędna jest modyfikacja filtru Kalmana. Niech τ będzie zbiorem punktów w czasie dla których zmienna yt jest nieobserwowana, to jest τ = {t : yt = NaN}. Dla t ∈ τ równania filtru Kalmana mają następującą postać αt|t−1 = Fαt−1|t−1, (1.17) Pt|t−1 = FPt−1|t−1F ′ +Q, (1.18) αt|t = αt|t−1, (1.19) Pt|t = Pt|t−1. (1.20) Estymacja komponentu trendu µ została wykonana w programie MatLab na podstawie własnej implementacji modelu. Funkcja wiarogodności została zmaksymalizowana za po- mocą algorytmu wyżarzania symulowanego, a uzyskana wartość funkcji log-wiarygodności równa była 161310,20. Oszacowana wartość wariancji składnika losowego ε równa jest σ̂2 ε =0,026798, przy odchyleniu standardowym tego oszacowania 4,0251·10−9. Bardzo zbliżone wartości oszacowań parametrów potwierdzają wniosek wyciągnięty w podroz- dziale 1.4 o istnieniu niestacjonarnego komponentu w danych o zapotrzebowaniu na moce przesyłowe krajowego systemu elektroenergetycznego. 2.1. MODEL HOLTA-WINTERSA 25 muje postać: ŷt+1|t = lt + bt + st−m, (2.2) lt = α(yt − st−m) + (1− α)(lt + bt), (2.3) bt = β(lt − lt−1) + (1− β)bt−1, (2.4) st = γ(yt − lt−1 − bt−1) + (1− γ)st−m, (2.5) gdzie l jest komponentem poziomu, b - komponentem trendu, s - komponentem sezo- nowym, natomiast m jest stałą równą okresowi sezonowości. Model szacowany jest za pomocą minimalizacji reszt min α,β,γ { T∑ i=1 ( yt − ŷt|t−1 )2 } = min α,β,γ { T∑ i=1 (yt − (lt−1 + bt−1 + st−m))2 } , przy ograniczeniu 0 < α, β, γ < 1. Przystępny opis modeli wygładzania wykładniczego znaleźć można w Hydman i Athanasopoulos (2014). Model postaci (2.2)-(2.5) oszacowany został w środowisku R przy użyciu pakietu forecast. Zakres danych pokrywał dwuletni okres 2011-2012, czyli 70176 kwadranse. Wybór podyk- towany był czasem estymacji. Oszacowania parametrów modelu i kryteria informacyjne przedstawione są w tablicy (2.1). Tablica 2.1: Oszacowania modelu Holta-Wintersa α β γ AIC BIC loglik 0,92 0,43 0,01 54152,34 54186,21 -27072,17 Źródło: Opracowanie własne. Reszty w liniowym gaussowskim modelu powinny charakteryzować się trzema właściwo- ściami. Commandeur i Koopman8 wymieniają je w malejącej skali ważności, tj. nie- zależność, homoskedastyczność i normalność. Sumaryczny opis testów diagnostycznych znajduje się w aneksie. Tablica 2.2: Testy diagnostyczne modelu Holta-Wintersa statystyka wartość p-value założenie niezależność Q(96) 24905,550 0 − r(1) 0,878 0,349 + r(96 ∗ 7) 24905,550 0 − homoskedastyczność H 1,034 0,005 − normalność N 7600,654 0 − Źródło: Opracowanie własne. 8Commandeur J.J.F., Koopman S.J., 2007, An Introduction to State Space Analysis, Oxford, str. 14. 26 ROZDZIAŁ 2. ANALIZA SEZONOWOŚCI Publikowanie diagnostyki reszt nie jest powszechne wśród autorów artykułów, w tym tych poświęconych modelowaniu popytu na energię. Niespełnienie najważniejszego założenia świadczy o występowaniu autokorelacji reszt, czyli nie wytłumaczonej przez model części wahań. Przyczyną może być m.in. występowanie w szeregu dodatkowej wersji sezono- wości. Rozwiązaniem problemu autokorelacji reszt może być uwzględnienie dodatkowych sezonowości, do czego potrzebny jest bardziej skomplikowany model. 2.2 Model Taylora Problem występowania dwóch rodzajów fluktuacji czasowych został dostrzeżony w lite- raturze przez innych autorów. W dziedzinie prognozowania zapotrzebowania na energię elektryczną tę tematykę poruszali Harvey i Koopman9, Harvey i inni10, Pedregal i Young11 oraz Dordonnat i inni12. Taylor13 przedstawił model będący rozszerzeniem klasycznego modelu Holta-Wintersa o drugą zmienną sezonową. Na podstawie danych półgodzinnych o popycie na energię elektryczną zauważył wyższość swojego modelu nad standardowym modelem Holta-Wintersa. Model zaproponowany przez Taylora, czyli model Holta-Wintersa z podwójną sezonowo- ścią, przyjmuje postać: ŷt+1|t = lt + bt + s (1) t−m1 + s (2) t−m2 , (2.6) lt = α(yt − s(1) t−m1 − s(2) t−m2 ) + (1− α)(lt + bt), (2.7) bt = β(lt − lt−1) + (1− β)bt−1, (2.8) s (1) t = γ(yt − lt−1 − bt−1 − s(2) t−m2 ) + (1− γ)s (1) t−m1 , (2.9) s (2) t = ω(yt − lt−1 − bt−1 − s(1) t−m1 ) + (1− ω)s (2) t−m2 . (2.10) Pakiet forecast do środowiska R umożliwia oszacowanie także modelu Taylora podwójnej sezonowości. Jego oszacowania przedstawione są w tablicy (2.3). 9Harvey A.C., Koopman S.J., 1993, Forecasting Hourly Electricity Demand Using Time-varying Spli- nes, Journal of the American Statistical Assocation, 88, 1228-1237. 10Harvey A.C., Koopman S.J., Riani M., 1997, The Modeling and Seasonal Adjustment of Weekly Observations, Journal of Business & Economic Statistics, Vol. 15(3), s. 354-368. 11Pedregal D.C., Young P.C., 2006, Modulated Cycles, an Approach to Modelling Periodic Components from Rapidly Sampled Data, International Journal of Forecasting, 22 (1), s. 181-194. 12Dordonnat V., Koopman S.J., Ooms M., Dessertaine A., 2008, An Hourly Periodic State Space Model for Modeling French National Electricity Load, 13Taylor J.W., 2003, Short-Term Electricity Demand Forecasting Using Double Seasonal Exponential Smoothing, Journal of Operational Research Society, 54, s. 799-805. 2.3. MODEL HYNDMANA I INNYCH 27 Tablica 2.3: Oszacowania modelu Taylora α β γ ω 0,704 2,35·10−9 0,848 0,679 Źródło: Opracowanie własne. Analiza reszt, zaprezentowana w tablicy 2.4, zaprzecza hipotezie o braku autokorelacji reszt. Oprócz występowania istotnej statystycznie autokorelacji pierwszego rzędu, reszty z modelu wykazują korelację ze swoimi opóźnionymi o tydzień wartościami. Rozszerzenie modelu o drugą sezonowość nie rozwiązało problemu autokorelacji 672-rzędu. Niespełnie- nie założenia o homoskedastyczności wraz z wnioskiem o autokorelacji oznacza, że należy spróbować znaleźć lepszą specyfikację modelu, dla której ten problem nie będzie wystę- pował. Tablica 2.4: Testy diagnostyczne modelu Taylora statystyka wartość p-value założenie niezależność Q(96) 69569,209 0 − r(1) 12271,640 0 − r(96 ∗ 7) 69569,209 0 − homoskedastyczność H 1,138 0 − normalność N 3582,789 0 − Źródło: Opracowanie własne. Jak zauważyli Hyndman i inni14, niektóre szeregi czasowe, jak na przykład liczba roz- mów przychodzących w bankowości detalicznej mierzona w 5-minutowych przedziałach czasu, może przejawiać dodatkowe roczne wahania sezonowe. Kolejnymi przytoczonymi tam przykładami są: dzienna liczba przyjęć szpitalnych, wypłaty gotówkowe z bankoma- tów, wykorzystanie energii elektrycznej czy zapotrzebowanie na wodę. Daje to podstawę, wraz z wynikami weryfikacji założeń, do analizowania występowania dodatkowego rodzaju wahań sezonowych. 2.3 Model Hyndmana i innych Taylor15 przedstawił propozycję uwzględniania trzech rodzajów sezonowości w modelowa- niu szeregów czasowych. Jakkolwiek model ten jest łatwo zrozumiały i względnie nieskom- plikowany w interpretacji, to wydaje się, że jego słabym punktem jest brak możliwości 14Hyndman R.J., De Livera A. M., Snyder R.D., 2010, Forecasting Time Series with Complex Seasonal Patterns Using Exponential Smoothing, Monasch University, Working Paper 15/09, str. 3. 15Taylor J.W., 2010, Triple Seasonal Methods for Short-term Load Forecasting, European Journal of Operational Research, 204, s. 139-152. 30 ROZDZIAŁ 2. ANALIZA SEZONOWOŚCI chastycznej sezonowości. Obserwowana zmienna jest sumą komponentów trendu, trzech komponentów sezonowych i błędu losowego. Oryginalnym pomysłem na modelowanie trzech rodzajów sezonowości w 15-minutowych danych popytu na energię elektryczną był model postaci: yt = µt + γt + θt + δt + εt, (2.19) µt = µt−1 + νt−1, (2.20) νt = νt + ζt, (2.21) γt = −γt−1 − γt−2 − . . .− γt−95 + ωt, (2.22) θt = φt (−θt−1 − θt−2 − . . .− θt−6) + (1− φt)θt + ηt, (2.23) δt = ϕt (−δt−1 − δt−2 − . . .− δt−11) + (1− ϕt)δt + ξt, (2.24) dla i = 1, 2, . . . , T , gdzie µ jest komponentem poziomu, ν - nachylenia, a obydwa mo- delują stochastyczny trend, natomiast komponenty γ, θ, δ reprezentują wahania sezonowe przypisane odpowiednio poszczególnym kwadransom, dniom i miesiącom. Błąd losowy ζt odpowiada za zmiany nachylenia komponentu trendu. Wzrostowi wariancji σ2 ζ odpowiada zwiększenie stochastycznych zmian w trendzie, natomiast w przypadku gdy σ2 ζ = 0 trend przyjmuje postać liniowego trendu deterministycznego. Błędy losowe ω, η, ξ, w równa- niach (2.19)-(2.24), są komponentami nieregularnymi, czyli białymi szumami o rozkładzie normalnym z zerowymi średnimi i wariancjami σ2 ω, σ 2 η, σ 2 ξ . Niestety, już w przypadku estymacji modelowania samej tylko sezonowości przypisanej do poszczególnych kwadransów, macierz przejścia w filtrze Kalmana ma rozmiar 95 × 95. Przy liczbie obserwacji 70176, dla danych za lata 2011-2012, wyznaczenie wartości funkcji wiarogodności dla jednego punktu trwa 4,5 minuty, przez co proces maksymalizacji jest zbyt długotrwały. Z tego powodu, zgodnie z poprzednim podrozdziałem, w części prezentowanych modeli zastosowana została harmoniczna forma sezonowości. Wówczas model (2.19)-(2.24) przyjmuje postać: yt = µt + k∑ i=1 γi,t + θt + δt + εt, (2.25) µt = µt−1 + νt−1, (2.26) νt = νt + ζt, (2.27) γi,t = γi,t−1 cos (λi) + γ∗i,t−1 sin (λi) + ωi,t, (2.28) γ∗i,t = −γi,t−1 sin (λi) + γ∗i,t−1 cos (λi) + ω∗i,t, (2.29) θt = φt (−θt−1 − θt−2 − . . .− θt−6) + (1− φt)θt + ηt, (2.30) δt = ϕt (−δt−1 − δt−2 − . . .− δt−11) + (1− ϕt)δt + ξt, (2.31) dla i = 1, 2, . . . , k, gdzie k jest liczbą harmonik niezbędną do przybliżenia fluktuacji 15-minutowych. Komponent γ∗ pełni tą samą rolę co s∗(i) w modelu (2.11)-(2.18). W 2.4. MODELE ZE STOCHASTYCZNĄ I DETERMINISTYCZNĄ SEZONOWOŚCIĄ31 praktyce, jak zauważył Harvey20, zakłada się, że Var(ωi,t) = Var(ω∗i,t) = σ2 i = σ2 ω, dla i = 1, 2, . . . , k, co pozwala istotnie zredukować czas estymacji przy niewielkiej stracie w dopasowaniu. 2.4.1 Szczegółowa specyfikacja modelu Zanim zaprezentuję wyniki modeli z różnymi rodzajami sezonowości chciałbym zwrócić uwagę na istotne szczegóły estymacji. W modelu (2.19)-(2.24) niezbędne jest powiązanie parametrów wariancji σ2 ω i σ2 ζ . Jeśli nie, to może się zdarzyć, że zwiększająca się wa- riancja ωt może pomniejszać wariancję błędu εt i tym samym przejmować część wahań błędu z równania obserwacyjnego. Innymi słowy, wariancja przypisywana komponentowi nieregularnemu, czyli błędowi losowemu z równania obserwacyjnego może dążyć do zera, σ2 ε −→ 0. Nie jest to pożądane, ponieważ nie można w takim przypadku określić jaka część wahań ma być przypisywana stochastycznemu komponentowi sezonowości a jaka błędowi losowemu w równaniu pomiarowym. W związku z tym, relacje parametrów qω = σ2 ω/σ 2 ε oraz qζ = σ2 ζ/σ 2 ε muszą być ustalone na odpowiednim poziomie. Relacje te mogą być wy- znaczone na podstawie funkcji wzmocnienia filtra dla komponentów sezonowych i trendu. Opracowania wykorzystujące modele z trendem i sezonowością nie zawierają wielu przy- kładów dyskusji nad relacją między błędami losowymi. Problem relacji między warian- cjami σ2 ω i σ2 ζ pojawia się w publikacjach dotyczących filtru Hodricka-Prescotta, tak jak na przykład w Harvey’u i Jaegerze21 czy Harvey’u i Trimburze22. W związku z tym zdecydowałem przedstawić jak ustalenie zależności między tymi parametrami wpływa na własności modelu oraz zaprezentować strategię estymacji modeli z trendem i sezonowością. Wyznaczenie wpływu relacji między wariancjami błędów losowych rozpocznę od rozwa- żenia modelu (2.19)-(2.24) ograniczonego do jednego rodzaju wahań sezonowych postaci: yt = µt + γt + εt, (2.32) µt = µt−1 + νt−1, (2.33) νt = νt + ζt, (2.34) γt = −γt−1 − γt−2 − . . .− γt−s+1 + ωi,t, (2.35) gdzie t = 1, 2, . . . , T . Dla uproszczenia wywodu załóżmy miesięczną częstotliwość obser- 20Harvey A. C., 1989, Forecasting, Structural Time Series Models and the Kalman Filter, Cambridge University Press, Cambridge, str. 126 oraz str. 180. 21Harvey A.C., Jaeger A., 1993, Detrending, Stylized Facts and the Business Cycle, Journal of Applied Econometrics, 8(3), s. 231-247. 22Harvey A.C., Trimbur T., 2008, Trend Estimation and the Hodrick-Prescott Filter, Journal of the Japan Statistical Society, 38(1), s. 41–49. 32 ROZDZIAŁ 2. ANALIZA SEZONOWOŚCI wowanej zmiennej, tj. s = 12. Wówczas, komponent trendu można przedstawić jako: µt = ζt (1− L)2 , a funkcję generującą kowariancję następująco: gµ(L) = 1 (1− L)2 1 (1− L−1)2 σ2 ζ = hµ(L)σ2 ζ . Komponent sezonowy może być przedstawiony za pomocą równania: γ = (1− L) (1− L12) ωt, a jego funkcja generująca kowariancję jako23: gγ(L) = (1− L) (1− L12) (1− L−1) (1− L−12) σ2 ω = hγ(L)σ2 ω. Funkcję wzmocnienia filtra wyznacza się podstawiając L = eiλ. Funkcja zysku z filtra dla komponentu sezonowego to: GL = gµ(λ) σ2 ε + gµ(λ) + gγ(λ) = qζhµ(λ) 1 + qζhµ(λ) + qωhγ(λ) = qζ 1 (2−2 cosλ)2 1 + qζ 1 (2−2 cosλ)2 + qω (2−2 cos 12λ) (2−2 cosλ) , (2.36) GS = gγ(λ) σ2 ε + gµ(λ) + gγ(λ) = qωhγ(λ) 1 + qζhµ(λ) + qωhγ(λ) = qω (2−2 cos 12λ) (2−2 cosλ) 1 + qζ 1 (2−2 cosλ)2 + qω (2−2 cos 12λ) (2−2 cosλ) , (2.37) gdzie qζ = σ2 ζ σ2 ε oraz qω = σ2 ω σ2 ε . Funkcje wzmocnienia filtra mają inny przebieg w zależności od parametrów qζ i qω. Wykres 2.1 przedstawia funkcje wzmocnienia filtra dla qζ ustalonego na poziomie qζ = 0.001, oraz dla różnych wartości qω. 23Na podstawie wzoru (1 + L+ L2 + . . .+ Ls−1) = (1− L)/(1− L12). 2.4. MODELE ZE STOCHASTYCZNĄ I DETERMINISTYCZNĄ SEZONOWOŚCIĄ35 Pokazane zostało jak ustalenie relacji między błędami w równaniu trendu i sezonowości przekłada się na włąściwości modelu, więc dalej zaprezentuję koncepcję estymacji mo- delu klasy (2.32)-(2.35). Strategia ta w literaturze określana jest jako wyprowadzenie (concentrate out) jednego z parametrów z funkcji wiarogodności. W pracach, w których autorzy korzystali z filtru Kalmana, pomija się to zagadnienie, a zgodnie z najlepszą wie- dzą autora nie istnieje opracowanie w polskiej literaturze poświęcone temu zagadnieniu. Wyprowadzenie parametru z funkcji wiarogodności pozwala uzyskiwać właściwe oszaco- wania parametrów bez konieczności ustalania relacji wariancji. Problem ten nie pojawiał się przy szacowaniu modelu przestrzeni stanów dla filtru Hodricka-Prescotta, ponieważ w modelu tym relacja wariancji błędów z równań pomiarowego i stanu była ustalona. Wstępem do opisu zagadnienia wyprowadzania parametru jest przedstawienie logarytmu funkcji wiarogodności modelu przestrzeni stanów, który dany jest wzorem25 logL = −NT 2 log 2π − 1 2 T∑ t=1 log |Ft| − 1 2 T∑ t=1 v ′ tF −1 t vt, (2.49) gdzie vt = yt − yt|t−1, dla t = 1, 2, . . . , T , a Ft jest macierzą kowariancji. Dla przypadku modelu z jedną zmienną funkcja wiarogodności ma postać26 logL = −T 2 log 2π − 1 2 T∑ t=1 log ft − 1 2 T∑ t=1 ν ′ tf −1 t νt. (2.50) Wyprowadzenie wariancji błędu pomiaru polega na ustaleniu pozostałych parametrów wa- riancji w proporcji do niego. Wtedy Var(εt) = σ2 ∗ht, dla np.: ht = 1, natomiast macierze postaci przestrzeni stanów modelu i macierze filtru Kalmana muszą zostać przemnożone27 przez σ2 ∗. Harvey28 zauważył, że taka reparametryzacja jest korzystna, ponieważ wykorzy- stuje liniowości do zredukowania wymiaru przestrzeni parametrów. Filtr Kalmana może być wówczas zastosowany niezależnie od wyprowadzonej zmiennej. Reparametryzacja modelu uprawnia do zapisania wektora szacowanych parametrów w postaci ψ = [ψ ′ ∗, σ 2 ∗] ′, gdzie ψ∗ jest wektorem zawierającym n− 1 parametrów i przedstawienia funkcji (2.50) w postaci logL = −T 2 log 2π − T 2 log σ2 ∗ − 1 2 T∑ t=1 log ft − 1 2σ2 ∗ T∑ t=1 v2 t ft . (2.51) Z pochodnej cząstkowej (2.51) otrzymuje się estymator wyprowadzonego parametru σ2 ∗. ∂ logL ∂σ2 ∗ = 0 ⇔ σ̃2 ∗(ψ∗) = 1 T T∑ t=1 v2 t ft (2.52) 25Stosuję oznaczenia przyjęte w opisie modeli przestrzeni stanów i filtru Kalmana w aneksie. 26Zastosowanie małych liter podkreślać ma przypadek jednej zmiennej. 27Przy zastosowaniu konwencji przyjętej w opisie modeli przestrzeni stanów i filtru Kalmana zapisać należy Var(ζ) = σ2 ∗Q, Pt+1|t = σ2 ∗P ∗ t+1|t i ft = σ2 ∗f ∗ t 28Harvey A.C., 1989, Forecasting, Structural Time Series Models ..., str. 126. 36 ROZDZIAŁ 2. ANALIZA SEZONOWOŚCI Funkcję wiarogodności z wyprowadzonym parametrem uzyskuje się z (2.51) przy podsta- wieniu (2.52) logL = −T 2 (log 2π + 1)− 1 2 T∑ t=1 log ft − T 2 log 1 T T∑ t=1 v2 t ft . (2.53) Zaprezentowane dalej modele zostały oszacowane przy wyprowadzeniu parametru warian- cji błędu losowego z równania pomiarowego σ2 ε z funkcją wiarogodności daną równaniem (2.53). 2.4.2 Modele z zero-jedynkową sezonowością Analogicznie jak w podrozdziałach 2.1-2.3, modelowanie danych 15-minutowych rozpo- częte zostało od najprostszego modelu, który był stopniowo rozbudowywany. Przedsta- wione w tym podrozdziale modele ekonometryczne zostały oszacowane przy wykorzysta- niu wyprowadzania parametru z funkcji wiarogodności. Zmienną obserwowaną w tych modelach jest logarytm 15-minutowego zapotrzebowania na moce przesyłowe KSE. Wyjściowy model stochastycznego trendu z sezonowością przyjmuje formę yt = µt + γt + εt, (2.54) µt = µt−1 + νt−1, (2.55) νt = νt + ζt, (2.56) γt = −γt−1 − γt−2 − . . .− γt−95 + ωt, (2.57) dla t = 1, 2, . . . , T , gdzie µ jest komponentem poziomu, ν - nachylenia, γ reprezentuje dzienne wahania sezonowe, natomiast ε odpowiada błędowi losowemu. Komponenty ζ,ω są zaburzeniami losowymi z zerową średnią i wariancjami odpowiednio σ2 ζ i σ2 ω. Został on oszacowany w wersji ze stochastyczną sezonowością, ponieważ wyniki analiz modeli z podrozdziałów 2.1-2.3 przesądziły o tym, że determinisyczna wersja sezonowości nie jest właściwa. Oszacowania modelu prezentuje tablica 2.6 w kolumnie (1). Model ten został porównany ze swoimi rozszerzonymi wersjami za pomocą współczynnika determinacji R2, który, zgodnie z Harvey’em29, ma postać R2 = 1− SSE / T∑ t=2 ( ∆yt −∆y )2 , gdzie SSE = f ∑T t=1 ṽ 2 t = (T−d)σ̃2, d jest liczbą komponentów w modelu, ∆y jest średnią z przyrostów zmiennej y, natomiast f̄ jest rozwiązaniem stanu ustalonego f = limt→∞ ft. Współczynnik determinacji oraz wartość wariancji błędu pomiarowego wskazują na prze- 29Harvey A. C., 1989, Forecasting, Structural Time Series Models and the Kalman Filter..., str. 268. 2.4. MODELE ZE STOCHASTYCZNĄ I DETERMINISTYCZNĄ SEZONOWOŚCIĄ37 Tablica 2.6: Wyniki dla modeli różnych sezonowości zero-jedynkowych (1) (2) (3) Model ze Model ze Model ze stochastyczną stochastycznymi stochastycznymi sezonowością sezonowościami: seznowościami: 15-minutową 15-minutową 15-minutową, i dzienną dzienną i miesięczną 105σ2 ε 1,9994 1,9932 1,9930 105σ2 ζ 1,3843 1,3711 1,3787 105σ2 ω 0,1326 0,1361 0,1345 105σ2 η 0,0147 0,0005 105σ2 δ 2,2278 10−3 logL 2,2307 2,2305 2,2304 R2 0,6512 0,6521 0,6522 Źródło: Opracowanie własne. wagę modelu ze stochastyczną sezonowością. Podobnie jak model Taylora30 jest rozszerzeniem modelu Holta-Wintersa, tak model po- staci yt = µt + γt + θt + εt, (2.58) µt = µt−1 + νt−1, (2.59) νt = νt + ζt, (2.60) γt = −γt−1 − γt−2 − . . .− γt−95 + ωt, (2.61) θt = φt (−θt−1 − θt−2 − . . .− θt−6) + (1− φt)θt + ηt, (2.62) jest rozszerzeniem (2.54)-(2.57). Równanie (2.58) zawiera drugą zmienną sezonową θt, która pozwala modelować wahania charakterystyczne dla poszczególnych dni. Struktura autoregresyjna θt pozwala przypisywać wartości poszczególnym dniom w modelu wysoko- częstotliwościowym. Pomysł na modelowanie drugiej sezonowości w ten sposób nie został wykorzystany dotąd w literaturze. Zaletą tego podejścia jest brak konieczności sztywnego odwoływania się do konkretnych opóźnień, dzięki czemu nie występuje problem zmian czasu omawiany na początku rozdziału 2. Oszacowania parametrów tego modelu znaj- dują się w tablicy 2.6 w kolumnie (2). Zastosowana została stochastyczna wersja dziennej sezonowości, a wybór ten został potwierdzony przez niezerowe oszacowanie parametru σ2 η. 30Taylor J.W., 2003, Short-term Electricity Demand Forecasting ... 40 ROZDZIAŁ 2. ANALIZA SEZONOWOŚCI Tablica 2.7: Wyniki dla modeli różnych sezonowości trygonometrycznych (1) (2) (3) Model z dwoma Model z trzema Model z trzema rodzajami rodzajami rodzajami stochastycznymi stochastycznymi stochastycznymi sezonowościami sezonowościami sezonowościami z trendem z trendem z trendem stochastycznym stochastycznym deterministycznym 105σ2 ε 2,6916 0,9692 1,3343 105σ2 ζ 3,6358 0,0087 0∗ 105σ2 ω 3,4673 122,8860 139,1378 10−3 logL 2,1351 2,2258 2,2183 R2 0,5314 0,8210 0,7926 ∗ oznacza, że parametr nie był szacowany, tylko zostało na niego nałożone ograniczenie. Źródło: Opracowanie własne. danych został potwierdzony wynikami modelu z trzema rodzajami sezonowości z determi- nistycznym trendem. Testowana testem ilorazu wiarygodności hipoteza H0 : σ2 ζ = 0 także musi zostać odrzucona na poziomie istotności 5%36. Stwierdzić zatem należy, że model ze potrójną stochastyczną strukturą sezonową i stochastycznym trendem jest właściwy do modelowania analizowanego szeregu. Przedstawione powyżej analizy potwierdzają założoną dodatkową hipotezę o istnieniu potrójnej struktury sezonowej w szeregu zapotrzebowania na moce przesyłowe krajowego systemu elektroenergetycznego. Dodatkowo pokazałem, że wahania sezonowe dzienne, tygodniowe i roczne powinny być modelowane w sposób stochastyczny tak jak komponent trendu. Rozdziały 1 i 2 zawierają analizy składowych długo- trendu i krótkookresowych, do których zalicza się sezonowość, natomiast pozostała część wahań, czyli komponent średniookresowy jest przedmiotem rozdziału 3. 36Wartość statystyki testowej była równa 1500, natomiast wartość krytyczna χ2 1 =3,84 Rozdział 3 Wpływ cyklu koniunkturalnego Zgodnie z teorią analizy szeregów czasowych, każdy z nich możemy zdekomponować na składniki długo-, średnio- oraz krótko-okresowe. Jak już zostało to zaprezentowane w pierwszym rozdziale, w szeregu zapotrzebowania na moce przesyłowe krajowej sieci elek- troenergetycznej występuje składowa trendu. Analizie komponentu sezonowego, zalicza- nego do wahań krótkookresowych, poświęcony został natomiast rozdział drugi. Składnik średniookresowy utożsamiany jest z tą częścią wahań, która związana jest z cyklem ko- niunkturalnym, a właściwe jego uwzględnienie i modelowanie może być przydatne również w wyjaśnianiu zmienności szeregu wysokoczęstotliwościowego i w prognozowaniu jego war- tości w średnim horyzoncie czasowym. Rysunek 3.1: Logarytm krajowego zużycia energii elektrycznej i jego prognoza Źródło: Opracowanie własne. Rozważmy miesięczny szereg czasowy krajowego zużycia energii elektrycznej dla okresu od stycznia 1993 do sierpnia 2013 roku. Dla próbki obserwacji 1993m01-2007m12, osza- 41 42 ROZDZIAŁ 3. WPŁYW CYKLU KONIUNKTURALNEGO cowany został model SARIMA1 dla logarytmu zmiennej, a następnie na jego podstawie wyznaczono prognozę na następne 24 miesięce. Rysunek 3.1 przedstawia wyznaczone pro- gnozy w porównaniu z rzeczywiście zaobserwowanymi wartościami. Prognozowana war- tość zużycia znacznie przewyższa zaobserwowaną. Widocznie mniejsze niż przewidywane zużycie może być związane z kryzysem zapoczątkowanym przez problemy na amerykań- skim rynku kredytów hipotecznych. Wydaje się, że wiedza o stanie gospodarki, o okresach prosperity i recesji, mogłaby istotnie poprawić prognozy popytu na energię elektryczną. W związku z tym, poniżej zaprezentowana została weryfikacja hipotezy, że przy wyko- rzystaniu informacji o wahaniach średniookresowych stanu gospodarki uzyskuje się lepsze prognozy popytu na energię elektryczną. Wyliczyłem nie tylko dwie zmienne ilustrujące stan koniunktury, ale także przedstawiłem propozycję szeregu, za pomocą którego pro- gnozować można stan gospodarki w przyszłości. Pokazałem w dalszej części rozdziału, że różne szeregi popytu na energię elektryczną zawierają składowe średniookresowe identy- fikowane z cyklem koniunktralny. Przedstawione poniżej próby i wyniki skwantyfikowa- nia polskiego cyklu koniunkturalnego zostały dalej wykorzystane do zbadania istnienia zależności przyczynowo-skutkowej z szeregami opisującymi zapotrzebowanie na energię elektryczną oraz jej ewentualnego kierunku. Ostatnia część rozdziału poświęcona jest modelowaniu i prognozowaniu danych wysokoczęstotliwościowych z wykorzystaniem wa- hań cyklicznych jako zmiennej egzogenicznej. 3.1 Równoczesny i wyprzedzający wskaźnik koniunk- tury Najpowszechniej przytaczana definicja wahań koniunkturalnych pochodzi od Mitchella2. Cykl koniunkturalny, zgodnie z nią, jest częścią wahań zagregowanej aktywności gospo- darczej w kraju, w którym praca zorganizowana jest w przedsiębiorstwach, przy czym uwzględnia się tylko te rodzaje aktywności, które mają systematyczny i ekonomiczny cha- rakter. Cykl składa się z okresu ekspansji przejawiającego się w tym samym czasie w wielu typach działalności gospodarczej, okresu dekoniunktury, recesji i ożywienia łączącego się z okresem ekspansji kolejnego cyklu. Sekwencja tych zmian ma charakter powtarzalny. Czas trwania cyklu waha się od ponad roku do dziesięciu, dwunastu lat. Poszczególne cykle nie mogą być dzielone na krótsze o podobnym charakterze. Według definicji Gau- dreault i innych3 termin „cykl koniunkturalny” odnosi się do wspólnych ruchów w bardzo 1Liczby opóźnień w modelu zostały wybrane za pomocą procedury od ogółu do szczegółu i miały postać odpowiednio SARIMA (3, 0, 0)× (0, 1, 1)12. Model został oszacowany w programie Stata. 2Mitchell W.C., 1927, Business Cycles: The Problem and Its Setting, National Bureau of Economic Research, Nowy Jork, s. 481. 3Gaudreault C., Lamy R., Liu Y., 2003, New Coincident, Leading and Recession Indexes for the Canadian Economy: An Application of the Stock and Watson Methodology, Department of Finance, Working Paper 2003-12, str. 2 3.1. RÓWNOCZESNY I WYPRZEDZAJĄCY WSKAŹNIK KONIUNKTURY 45 ponent ψn,t jest cyklem stochastycznym rzędu n, dla dodatnich n, jeśli[ ψi,t ψ∗i,t ] = ρ [ cosλc sinλc − sinλc cosλc ][ ψi,t−1 ψ∗i,t−1 ] + [ ψi−1,t−1 ψ∗i−1,t−1 ] , (3.4)[ ψ1,t ψ∗1,t ] = ρ [ cosλc sinλc − sinλc cosλc ][ ψ1,t−1 ψ∗1,t−1 ] + [ κt κ∗t ] , (3.5) gdzie i = 2, . . . , n, κt ∼ N(0, σ2 κ), κ∗t ∼ N(0, σ2 κ∗). Parametr ρ jest parametrem tłumienia (damping parameter) i spełnia 0 < ρ ≤ 1. Model wykorzystany do estymacji wskaźników koniunktury jest wielorównaniowym mo- delem strukturalnego szeregu czasowego ze wspólnym stochastycznym cyklem. Reprezen- tacja modelu w postaci przestrzeni stanów składa się z dwóch równań, pierwszego, obser- wacyjnego dla yt, drugiego dla zmiennych stanu opisującego ich dynamikę. Wektor zmien- nych stanu αt = (µ′t, ψ ′ t) ′, gdzie µ′t = (µ′m,t, µ ′ m−1,t, . . . , µ ′ 1,t) i ψ′t = (ψ′n,t, ψ ′ n−1,t, . . . , ψ ′ 1,t), gdzie µ′j,t = (µ1′ j,t, µ 2′ j,t, . . . , µ N ′ j,t ), gdzie superskrypt oznacza komponent trendu j-tego rzędu dla n-tej obserwowalnej zmiennej. Równanie pomiarowe ma postać yt = z′tαt + εt, dla t = 1, . . . , T , gdzie pierwszymi (m + 1) elementami wektora zt są jedynki, a zera pozostałymi. Równanie przejścia (stanu) dla komponentów trendu jest postaci µt = Umµt−1 + imζt, gdzie Um jestm×m wymiarową, górnotrójkątną macierzą jedynek, natomiast im, jestm×1 wektorem jedynek. Macierz kowariancji wektora zaburzeń jest równa V ar(itζt) = imi ′ mσ 2 ζ . Równanie przejścia dla komponentów cyklicznych jest następujące ψt = Tψψt−1 + in ⊗ [ κt κ∗t ] , gdzie macierz przejścia Tψ może być zapisana w postaci Tψ = In ⊗ ρ [ cosλc sinλc − sinλc cosλc ] + Sn ⊗ I2, gdzie Sn jest n × n wymiarowa macierzą zer z wyjątkiem elementów diagonalą równym jedynkom, a I2 jest macierzą identycznościową wymiaru 2×2, natomiast⊗ oznacza iloczyn Kroneckera. Bardziej szczegółowy opis modelu i jego zapisu w postaci przestrzeni stanów znajduje się w Harvey i Jaeger11, Harvey i Trimbur12 oraz Azevedo i inni13. 11Harvey A.C., Jaeger A., 1993, Detrending, Stylized Facts... 12Harvey A.C., Trimbur T.M., 2003, General Model-Based Filters... 13Azevedo J.V.e, Koopman S.J., Rua A., 2003, Tracking Growth and the Business Cycle: A Stochastic 46 ROZDZIAŁ 3. WPŁYW CYKLU KONIUNKTURALNEGO Po przedstawieniu strukturalnych modeli szeregów czasowych przedstawię parametry mo- delu tego typu, który wykorzystałem do estymacji równoczesnego wskaźnika koniunktury. Wybór m = 2 i k = 6 ten został uzasadniony przez Azevedo i inni14 tym, że wówczas model ten ma właściwości zbliżone do filtru Baxtera-Kinga. Model ten równoległego wyestymowany był na N = 5 zmiennych, a parametry były estymowane za pomocą me- tody największej wiarogodności. Częstotliwość λc wybrana została na poziomie π/48, który odpowiada wahaniom 96-miesięcznym i jest uzasadniony w literaturze przez Aze- vedo i innych15 jako odpowiednia wartość do modelowania cyklu koniunkturalnego w Unii Europejskiej i Stanach Zjednoczonych. W modelu wskaźnika równoczesnego ogranicze- nie na parametr ładunku przy szeregu produkcji premysłowej zostało narzucone tak, że δ1 = 1. W modelu indeksu wyprzedzającego podobne ograniczenie zostało nałożone na komponent cykliczny szeregu różnicy między rentownością 5-letnich obligacji skarbowych a 3-miesięcznym WIBORem. Model stochastycznego trendu (3.2)-(3.3) może być postrze- gany jako filtr dolnoprzepustowy, czyli pozostawiający z całego zakresu wahań tylko te o najniższych częstotliwościach, czyli takie, których okres jest bardzo duży lub nieskoń- czony. Model (3.4)-(3.5) jest filtrem średnioprzepustowym, gdyż z obserwowanego szeregu pozostawia tylko średnioczęstotliwościowe oscylacje. Harvey i Trimbur16 określają model (3.1)-(3.4) jako uogólniony dolno- i średnioprzepustowy filtr Butterwortha rzędu (m,n). Zgodnie z Harvey i Trimbur17, położenie maksimum funkcji wzmocnienia filtra średnio- przepustowego jest regulowane przez częstotliwość λ, natomiast kształt zależy od ustalenia parametru n. W oszacowanym modelu częstotliwość λ została ustalona na poziomie π/48, co odpowiada wahaniom 96-miesięcznym. Celem niniejszego podrozdziału nie jest jednak dokładne prezentowanie teorii filtracji, tylko zaprezentowanie modelu wynikającego z po- łączenia już istniejących. Tablica 3.2: Oszacowania modelu wskaźnika równoczesnego Parametry Zmienne Produkcja Handel przemysłowa detaliczny Import Eksport BWUK δi 1 0,434 1,742 1,728 0,733 σ2 i 0,016 0,033 0,297 0,300 0,077 δi/σ∆yi 43,434 14,597 42,092 46,637 13,496 Źródło: Opracowanie własne. Parametr ρ został oszacowany na poziomie ρ = 0.487. Tablica 3.2 zawiera oszacowania modelu wskaźnika równoczesnego., natomiast komponent cykliczny przedstawiony jest na Common Cycle Model for the Euro Area, Tinbergen Institute Discussion Paper, 2003-069/4. 14Azevedo J.V.e, Koopman S.J., Rua A., 2003, Tracking Growth..., str. 9. 15Azevedo J.V.e, Koopman S.J., Rua A., 2003, Tracking Growth..., str. 9. 16Harvey A.C., Trimbur T.M., 2003, General Model-Based Filters... 17Harvey A.C., Trimbur T.M., 2003, General Model-Based Filters..., str. 246 3.1. RÓWNOCZESNY I WYPRZEDZAJĄCY WSKAŹNIK KONIUNKTURY 47 rysunku 3.2. Zacienione obszary odpowiadają okresom od szczytu do dna koniunktury. Rysunek 3.2: Komponent cykliczny równoczesnego wskaźnika koniunktury 1990 1995 2000 2005 2010 2015 −6 −4 −2 0 2 4 6 ·10−2 Źródło: Opracowanie własne. Graficzne porównania składnika cyklicznego z wynikami z filtru Hodricka-Prescotta, zmo- dyfikowanego filtru Hodricka-Prescotta i Christiano-Fitzgeralda, przedstawione w Skrzyp- czyńskim18 i Skrzypczyńskim19, poświadczają duże podobieństwo przebiegu. W okresie, który może być porównany20 stwierdzam podobieństwo za pomocą przypadających w tych samych okresach szczytów koniunktury w latach 1998, 2000 i 2004 między przedstawio- nymi w Skrzypczyńskim21 a uzyskanym przeze mnie komponentem cyklicznym wskaźnika koniunktury. Wysoką zbieżność można zauważyć przy porównaniu przebiegu okresu oży- wienia i ekspansji w latach 1999-2000, gdzie charakterystyczne dodatkowe lokalne mak- simum jest obserwowane zarówno w przypadku wyników Skrzypczyńskiego22 jak i mo- ich. Adamowicz i inni23 zaprezentowali szeroką analizę synchronizacji polskiego cyklu 18Skrzypczyński P., 2008, Wahania aktywności gospodarczej w Polsce i strefie euro, Materiały i Studia nr 227, Narodowy Bank Polski. 19Skrzypczyński P., 2010, Metody spektralne w analizie cyklu koniunkturalnego gospodarki Polskiej, Materiały i Studia nr 252, Narodowy Bank Polski. 20W pracach Skrzypczyńskiego przedmiotem analiz były dane za lata 1995-2007. 21Skrzypczyński P., 2010, Metody spektralne w analizie cyklu..., str. 142, 143, 146, 147 i 148. 22Skrzypczyński P., 2010, Metody spektralne w analizie cyklu..., str. 148. Dodatkowe lokalne maksimum widoczne jest dla komponentów cyklicznych uzyskanych za pomocą filtrów Hodricka-Prescotta, modeli UCARIMA i SVAR. 23Adamowicz E., Dudek S., Pachnicki D., Walczyk K., 2009, Synchronizacja cyklu koniunkturalnego polskiej gospodarki z krajami strefy euro w kontekście struktury tych gospodarek w Raport na temat pełnego uczestnictwa Rzeczypospolitej Polskiej w trzecim etapie Unii Gospodarczej i Walutowej, Narodowy Bank Polski, Warszawa. 50 ROZDZIAŁ 3. WPŁYW CYKLU KONIUNKTURALNEGO regów o różnych częstotliwościach. Pierwszy jest opisany w Harvey35 i można go rozumieć jako interpolację zmiennej kwartalnej pryz uzciu ymiennej kumulujcej. Drugi sposób, opi- sany w Mariano i Murasawa36 jest rozszerzeniem modelu Stocka i Watsona37 o zmienną kwartalną. Dodatkowe rozszerzenie modelu Mariano i Murasawy o zmienną kumulującą zaproponował Nunes38. Moja propozycja jest rozszerzeniem istniejącego modelu wspólnego stochastycznego cy- klu, wykorzystanego do wyliczenia wskaźnika koniunktury w poprzednim podrozdziale, o zmienną kumulującą. Rozważmy przedstawiony w poprzednim podrozdziale model wspól- nego stochastycznego cyklu z pewnymi rozszerzeniami. Niech, podobnie jak poprzednio, yit będzie i-tą obserwowaną zmienną, która modelowana jest następująco: yit = µit + δi { cos (ξiλ)ψn,t + sin (ξiλ)ψ∗n,t } + εit, (3.6) µit = µit−1 + νit−1, (3.7) νit = νit−1 + ζit, (3.8)[ ψj,t ψ∗j,t ] = ρ [ cosλ sinλ − sinλ cosλ ][ ψj,t−1 ψ∗j,t−1 ] + [ ψj−1,t−1 ψ∗j−1,t−1 ] , (3.9)[ ψ1,t ψ∗1,t ] = ρ [ cosλ sinλ − sinλ cosλ ][ ψ1,t−1 ψ∗1,t−1 ] + [ κt κ∗t ] , (3.10) dla j = 2, . . . , n, gdzie t jest indeksem czasu t = 1, 2, . . . , T , natomiast i indeksem dla zmiennych i = 1, 2, . . . , N . Zgodnie z Azevedo i innymi39, parametr ξi interpretuje się jako przesunięcie fazowe i-tej zmiennej względem cyklu referencyjnego, przy czym ξi > 0 (ξi < 0) oznacza wyprzedzanie (opóźnianie) względem cyklu szeregu referencyjnego. Za- gadnienie dodania do modelu zmiennej o niższej częstotliwości wiąże się z dodaniem kolej- nej zmiennej stanu, która jest nieobserwowalnym odpowiednikiem tej zmiennej. Zmienna kwartalna jest obserwowana tylko w co trzecim miesiącu, więc niezbędne jest wprowadze- nie zmian do modelu W efekcie braku informacji o wartościach zmiennej kwartalnej w niektórych okresach czasu niezbędne jest wprowadzenie zmian do modelu (3.6)-(3.10). Niech yq będzie zmienną kwartalną, która ma zostać włączona do miesięcznego modelu określonego równaniami (3.6)-(3.10). Zmienna kwartalna określona jest jedynie w punk- tach czasu t = 3s, dla s = 1, 2, . . . , bT 3 c. Niech y†q będzie zmienną zdefiniowaną dla t = 1, 2, . . . , T tak, że y†q,3s = yq,3s i nieobserwowaną w punktach czasu 3s − 2 i 3s − 1. Zmienna y†q,t nie jest sprzeczna z modelem (3.6)-(3.10), gdyż jest określona dla tych samych 35Harvey A. C., 1989, Forecasting, Structural Time Series... 36Mariano R.S., Murasawa Y., A New Coincident Index of Business Cycles Based on Monthly and Quarterly Series, PIER Working Paper 02-014, Philadelphia 2002. 37Stock J.H., Watson M.W., A Probability Model of The Coincident Economic Indicators, Opubliko- wane w: G. Moore and K. Lahiri, editors. The Leading Economic Indicators: New Approaches and Forecasting Records, Cambridge 1990. 38Nunes L.C., 2005, Nowcasting Quarterly GDP Growth in a Monthly Coincident Indicator Model, Journal of Forecasting, 24, s. 575–592. 39Azevedo J.V.e, Koopman S.J., Rua A., 2003, Tracking Growth..., str. 5. 3.2. MIESIĘCZNE PKB 51 punktów w czasie, t = 1, 2, . . . , T . Równanie pomiarowe dla zmiennej y†q ma postać: yq,δt = µq,t + δqψn,t + εq,t. (3.11) Równanie (3.11) prawdziwe jest tylko dla punktów w czasie, w których dostępne są obser- wacje zmiennej kwartalnej, to jest w ostatnim miesiącu każdego kwartału. W przypadku nałożenia ograniczenia δq = 1, wspólny komponent cykliczny będzie miał amplitudę równą amplitudzie składnika cyklicznego zmiennej kwartalnej. Konsekwentnie, parametry δi przy zmiennych miesięcznych będą estymowane. Nieobserwowany odpowiednik zmiennej kwartalnej wyliczany jest ze wzoru y†t = µqm,t + ψn,t + εt, (3.12) czyli jako suma komponentu trendu, składnika cyklicznego i składnika nieregularnego. Metoda pozwala na interpolację brakujących wartości składnika cyklicznego zmiennej kwartalnej na podstawie zmiennych miesięcznych. Składnik cykliczny jest szeregiem od- chyleń od trendu w określonym zakresie częstotliwości. W analizach koniunktury na pod- stawie takiego wskaźnika możliwe jest wyznaczenie punktów zwrotnych cyklu, w których następuje zmiana średniookresowych tendencji. Punkty zwrotne, nazywane też punktami zmiany tendencji, są lokalnymi ekstremami szeregu reprezentującego koniunkturę. Na wy- znaczenie punktów zwrotnych nie ma wpływu monotoniczna transformacja szeregu. Fakt ten jest wykorzystany w Azevedo i inni40, gdzie składnik cykliczny został potraktowany jako zmienna stanu, przez co analiza koniunktury w strefie euro została dokonana na in- terpolacji składnika cyklicznego PKB. Wyliczenie miesięcznego PKB wymaga dodania do modelu ograniczenia sumowalności wartości miesięcznych do sum kwartalnych za pomocą zmiennej kumulującej, którą prezentuję poniżej. Niech ypkb,3s będzie kwartalnym szeregiem produktu krajowego brutto obserwowanym tylko dla miesięcy kończących kwartał. Zgodnie z powyższym zdefiniowana musi zostać zmienna y†pkb,t, która jest nieobserwowanym odpowiednikiem ypkb,3s dla wszystkich okresów czasu t = 1, 2, . . . , T . Zagadnienie miesięcznego PKB polega natomiast na wyliczeniu y†pkb,t, takiego, że ypkb,3s = y†pkb,3s−2 + y†pkb,3s−1 + y†pkb,3s. (3.13) Rozwiązaniem problemu sumowalności miesięcznego PKB do wartości kwartalnych jest użycie zmiennej kumulującej, tak jak przedstawiono to w Harvey’u41. Zmienna kumu- lująca, yft , dla zmiennej przepływu określona jest w sposób, że dla pierwszego miesiąca kwartału równa jest nieobserwowanemu miesięcznemu PKB, y†pkb,t, natomiast dla kolej- nych miesięcy równa jest sumie yft−1, swojej przeszłej wartości, i y † pkb,t, czyli miesięcznego PKB w tym okresie. Rekurencyjna zależność dla zmiennej kumulującej może być zapisana 40Azevedo J.V.e, Koopman S.J., Rua A., 2003, Tracking Growth... 41Harvey A. C., 1989, Forecasting, Structural Time Series..., str. 313. 52 ROZDZIAŁ 3. WPŁYW CYKLU KONIUNKTURALNEGO w postaci yf3(s−1)+r = ϕ3(s−1)+ry f 3(s−1)+r−1 + y†3(s−1)+r, (3.14) gdzie r = 1, 2, 3, natomiast ϕ3(s−1)+r zdefiniowane jest jako ϕ3(s−1)+r = { 0 dla r = 1, 1 dla r = 2, 3. (3.15) W przypadku danych kwartalnych w miesięcznym modelu ϕδ(t−1)+s jest indykatorem, który przyjmuje wartość zero dla pierwszego miesiąca w kwartale. Zastosowany model miesięcznego PKB jest oparty na równaniach (3.6)-(3.10), (3.13)- (3.14). Wektor zmiennych obserwowanych ma postać y =  Produkt krajowy brutto Produkcja przemysłowa Import Eksport Handel detaliczny WIG20 Rentowność 5-letnich obligacji skarb. – 3-mies. WIBOR Bieżący wskaźnik ufności konsumenckiej  . (3.16) Ze względu na identyfikację parametrów musi zostać nałożone ograniczenie δpkb = 1, przez co wspólny komponent cykliczny będzie miał amplitudę równą amplitudzie składnika cyklicznego zmiennej kwartalnej. Zastosowany został model stochastycznego cyklu rzędu 6 oraz model lokalnego liniowego trendu. Wybór λ = π 48 , m = 2 oraz n = 6 uzasadniony jest analogicznie jak w podrozdiale 3.1. W modelu rozszacowywany był odsezonowany szereg produktu krajowego brutto w cenach stałych przedstawiony jako indeks o bazie w roku 1995 publikowany zgodnie z metodologią ESA95. Oszcowania parametrów modelu znajdują się w tablicy 3.4. Komponent cykliczny miesięcznego PKB ma bardzo zbliżony przebieg do komponentu cyklicznego z modelu wskaźnika równoczesnego. Wnioski wyciągnięte z porównania wa- hań średniookresowych wskaźnika koniunktury pozostają niezmienione dla komponentu średniookresowego miesięcznego PKB. Zgodnie z terminologią NBER punktem zwrotnym lub punktem zmiany tendencji (tur- ning point) jest pojedynczy okres czasu, w którym następuje zmiana pomiędzy okresami wzrostu i kurczenia się gospodarki. Pikiem lub szczytem koniunktury (peak) jest punkt czasu, w którym gospodarka przestaje się rozwijać i zaczyna się kurczyć. Odwrotnym przypadkiem jest dno (trough), które jest okresem granicznym między fazą kurczenia się a fazą wzrostu. Zgodnie z definicją cyklu koniunkturalnego okres ekspansji kończy się szczytem i po nim zaczyna się okres dekoniunktury, czyli kurczenia się zagregowanej 3.3. DEKOMPOZYCJA SZEREGÓW ENERGII ELEKTRYCZNEJ 55 Tablica 3.5: Porównanie punktów zwrotnych wskaźników i PKB MPKB Wskaźniki OECD Kwartalne COIN LEI Opublikowane LEI PKB 1998.04 S -3 0 -1 -9 1998q2 1998.12 D +1 -1 - -8 1999q1 2000.09 S +1 -6 - -18 2000q1 2002.04 D -1 +11 +8 -10 2002q4 2004.04 S 0 -1 0 -1 2004q2 2005.04 D 0 0 +2 0 2005q2 2008.02 S 0 -6 0 -12 2008q2 2009.05 D -1 -3 +3 -7 2010q1 2011.04 S +8 0 +7 -8 2011q3 2013.02 D 0 -8 - -3 - W tabeli przyjęto oznaczenia: S - szczyt, D - dno koniunktury. MPKB - oznacza mie- sięczne PKB, COIN - wskaźnik koniunktury, LEI - indeks wyprzedzający. Źródło: Opra- cowanie własne. zgodność oszacowań. Zauważalne na rysunkach 3.3 i 3.4 krótsze fazy spowolnienia i rece- sji niż ożywienia i rozkwitu są potwierdzone przez analizę punktów zwrotnych, co zostało dostrzeżone także w Gradzewiczu i innych47. Do formalnego porównania punktów zwrot- nych wykorzystałem statystykę CCcorr przedstawioną przez Artisa i innych48. Wartość uzyskanej statystyki równa 80,13% może być zinterpretowana jako istnienie silnej korelacji między komponentem cyklicznym równoczesnego wskaźnika koniunktury a komponentem cyklicznym kwartalnego PKB49. 3.3 Dekompozycja szeregów energii elektrycznej Jednym ze sposobów weryfikacji hipotezy o wpływie cyklu koniunkturalnego na popyt na energię elektryczną może być wyodrębnienie komponentów cyklicznych z danych go opisu- jących i porównanie ich z cyklem koniunkturalnym. Składowa średniookresowa wyznaczy- łem za pomocą filtru średnio-przepustowego Christiano-Fitzgeralda, który był stosowany w analizach polskiego cyklu koniunkturalnego, na przykład przez Skrzypczyńskiego50 czy Gradzewicza i innych51. Drugim narzędziem wykorzystanym do sprawdzenia czy szeregi energii elektrycznej zawierają wahania w paśmie przypisywanemu wahaniom koniunktu- 47Gradzewicz M., Growiec J., Hagemejer J., Popowski P., 2010, Cykl koniunkturalny w Polsce – wnioski z analizy spektralnej, Bank i Kredyt, 41(5), s. 41-76. 48Artis M.J., Kontolenis Z.G., Osborn D.R., 1997, Business Cycles for G7 and European Countries, Journal of Business, 70(2), s. 249-279. 49Dla porównania komponentu cyklicznego miesięcznego PKB z komponentem cyklicznym kwartalnego PKB uzyskałem tą samą wartość statystyki CCcorr. 50Skrzypczyński P., 2010, Metody spektralne... 51Gradzewicz M., Growiec J., Hagemejer J., Popowski P., 2010, Cykl koniunkturalny... 56 ROZDZIAŁ 3. WPŁYW CYKLU KONIUNKTURALNEGO ralnym była analiza spektralna. Szeregi czasowe o popycie na energię elektryczną dla Polski zawarte są w seriach publikacji Agencji Rynku Energii52. Tabele z danymi w nich zawartymi, nie są dostępne w wersji elektronicznej, więc niezbędne było ich własnoręczne przepisanie. Powstała w ten sposób baza danych zawierała 2 szeregi miesięczne, 9 rocznych i 4 z łamaną częstotliwością. Do analizy tych ostatnich niezbędne było niestandardowe podejście. 3.3.1 Dane roczne Częstotliwość roczną mają szeregi całkowitego zużycia energii elektrycznej w kraju oraz w podziale na gałęzie gospodarki: przemysł, budownictwo, transport, pozostali odbiorcy. Dodatkowo analizowane były szeregi zużycia energii w lokalach mieszkalnych oraz w go- spodarstwach domowych i rolnych w podziale na miasto i wieś. Dane pokrywają okres od 1990 do 2012 roku, czyli dostępne były 23 obserwacje. Rysunek 3.5: Znormalizowane funkcje gęstości widmowej szeregów rocznych Źródło: Opracowanie własne. Rysunek 3.5 przedstawia znormalizowane funkcje gęstości widmowej dla szeregów rocz- nych. Zmienne zostały pozbawione liniowego trendu oraz przemnożone przez okno Blackmana- 52ARE, Informacja statystyczna o energii elektrycznej, Agencja Rynku Energii, Warszawa. ARE, Sta- tystyka Elektroenergetyki Polskiej, Agencja Rynku Energii, Warszawa. ARE, Sytuacja energetyczna w Polsce - krajowy bilans energii, Agencja Rynku Energii, Warszawa. ARE, Sytuacja w elektroenergetyce, Agencja Rynku Energii, Warszawa. 3.3. DEKOMPOZYCJA SZEREGÓW ENERGII ELEKTRYCZNEJ 57 Harrisa, co miało na celu rozwiązanie problemu wycieku mocy53. Maksima funkcji gęstości dla zużycia ogółem i w przemyśle przypadają na wahania powtarzające się odpowiednio co około 8 i 6 lat. Okresy wpadają w definicyjny zakres wahań koniunkturalnych, co może być argumentem potwierdzającym hipotezę o obecności wahań koniunkturalnych w szeregach popytu na energię elektryczną. W zużyciu energii elektrycznej przez podmioty, których działalność gospodarcza jest klasyfikowana jako transport wykazują silne wahania powtarzające się co 2,2 roku. Dla pozostałych szeregów dominujące fluktuacje mają okres zbliżony do górnej granicy definicji cyklu koniunkturalnego. Tablica 3.6: Wyniki analizy korelacji krzyżowych komponentów cyklicznych Zmienna Współczynnik Opóźnienie korelacji Zużycie w kraju 0,798 0 Przemysł 0,308 0 Budownictwo 0,212 0 Rolnictwo 0,302 -1 Transport 0,443 0 Pozostali odbiorcy 0,355 0 Lokale mieszkalne -0,355 0 Gosp. domowe i rolne na wsi -0,270 0 Gosp. domowe i rolne w mieście -0,202 1 Dodatnia wartość opóźnienia oznacza, że komponent cykliczny PKB jest opóźniony o daną liczbę okresów w stosunku do analizowanego szeregu. Źródło: Opracowanie własne. Przy użyciu filtra Christiano-Fitzgeralda wyznaczone zostały składowe średniookresowe rocznych szeregów energii elektrycznej i PKB. Tablica 3.6 przedstawia wyniki badania korelacji krzyżowej54 komponentów cyklicznych szeregów z cyklem produktu krajowego brutto. Odnotować należy silną korelację zużycia krajowego z PKB oraz umiarkowaną dodatnią korelację przemysłu a PKB. Komponenty cykliczne wszystkich zmiennych poza rolnictwem i miejskich gospodarstw domowych są równoczesne ze składową cykliczną PKB. Interesującym wnioskiem jest antycykliczność zapotrzebowania na energię elek- tryczną w lokalach mieszkalnych. Dane roczne w ogólności wykazują wysoką korelację z szeregiem produktu krajowego brutto, dlatego w następnym kroku sprawdzę istnienie korelacji dla dancyh miesięcznych. 3.3.2 Dane miesięczne Analogiczne postępowanie zostało przeprowadzone dla danych miesięcznych. W tym pod- rozdziale wykorzystałem szeregi produkcji energii elektrycznej w kraju oraz zużycia kra- 53Problem wycieku mocy został omówiony w podrozdziale 1.2. 54Korelacja krzyżowa jest narzędziem statystycznym służącym do analizowania korelacji między dwoma szeregami czasowymi z uwzględnieniem wzajemnych przesunięć w czasie. 60 ROZDZIAŁ 3. WPŁYW CYKLU KONIUNKTURALNEGO sięcznych w okresie 2003-2012, wobec czego wnioski powinny być prawidłowe. Lokalne maksima funkcji mocy przypadają dla częstości 2, 5, 13, 17, które odpowiadają waha- niom o okresach około 10, 4 lat, 18- oraz 14-miesięcznym. Pik dla wartości 20 odpowiada miesięcznej sezonowości. Wahania o częstotliwości 10-letniej są pozostałością po usunię- tym trendzie. Obecność wahań powtarzających się co 4 lata i co 1,5 roku świadczy o występowaniu wahań zgodnych z cyklem koniunkturalnym. Rysunek 3.7: Funkcje gęstości widmowej szeregów energii elektrycznej dostarczonej od- biorcom (a) ogółem, (b) na WN, (c) na SN i (d) na nN Źródło: Opracowanie własne. 3.4. KORELACJA I PRZYCZYNOWOŚĆ 61 3.4 Korelacja i przyczynowość W poprzednich podrozdziałach pokazane jest istnienie wahań o częstościach cyklu ko- niunkturalnego w wielu szeregach dotyczących energii elektrycznej. Niniejszy podroz- dział natomiast ma pokazać istnienie zleżności przyczynowo-skutkowej między energią elektryczną a stanem gospodarki. Najlepszym wyznacznikiem zagregowanego stanu go- spodarki jest produkt krajowy brutto. W związku z jego kwartalną częstotliwością zapro- ponowany został jego miesięczny odpowiednik. Miesięczna częstotliwość jest najwyższą, na jaką można bez większych zastrzeżeń rozszacować kwartalne PKB. Wyliczone w pod- rozdziałach 3.3.1 i 3.3.2 współczynniki korelacji pozwoliły badać istnienie przyczynowości. Średniookresowe komponenty produkcji i zużycia energii elektrycznej zostały wyodręb- nione z szeregów za pomocą filtru Christiano-Fitzgeralda. Wartości korelacji krzyżowej między tymi komponentami a składnikiem cyklicznym miesięcznego PKB jest najwyższa co do modułu dla 3 opóźnienia produkcji energii ogółem oraz 2 opóźnienia krajowego zużycia energii elektrycznej. Szeregi te wyprzedzają koniunkturę, co jest potwierdzane przez analizę punktów zwrotnych miesięcznych szeregów energii elektrycznej. Tablica 3.8: Punkty zwrotne szeregów produkcji i zużycia krajowego energii elektrycznej Szczyt/Dno Produkcja Zużycie MPKB D 1993.05 - - S 1996.02 1996.02 - D 1998.01 1996.11 1998.12 S 2000.09 1998.11 2000.09 D 2002.02 2002.02 2002.04 S 2004.04 - 2004.04 D 2005.01 - 2005.04 S 2006.01 2006.01 - D 2007.01 2006.12 - S 2007.11 2007.11 2008.02 D 2009.04 2009.04 2009.05 S 2011.05 2012.02 2011.04 D 2012.10 - 2013.02 Źródło: Opracowanie własne. Punkty zwrotne produkcji i zużycia energii elektrycznej, zaprezentowane w tablicy 3.8, potwierdzają wyniki korelacji krzyżowej. Dodatkowo, średnia długość cyklu w okresie 2002.01-2013.10 wynosi około 32 miesiące dla produkcji i 40 dla zużycia. Wartości te są zbliżone do długości cykli miesięcznego i kwartalnego PKB - odpowiednio 37 miesiące i 3 lata i 1 4 kwartału. Wartości te nie odbiegają bardzo od 2.35 roku uzyskanego z analizy widma danych 15-minutowych. Według definicji zaproponowanej przez Grangera60 stacjonarny proces xt jest przyczyną 60Granger C.W.J., 1969, Investigating Causal Relations by Econometric Models and Cross-Spectral 62 ROZDZIAŁ 3. WPŁYW CYKLU KONIUNKTURALNEGO innego procesu stacjonarnego yt, gdy xt jest lepiej prognozowany z wykorzystaniem yt niż bez niego. Zgodnie z definicją przedstawioną w Hamiltonie61 zmienna y nie jest przyczyną w sensie Grangera zmiennej x, jeśli średni błąd prognozy wartości x nie maleje przy wykorzystaniu opóźnionych wartości zmiennej y, innymi słowy gdy MSE [E (xt+s|xt, xt−1, . . . , yt, yt−1, . . .)] = MSE [E (xt+s|xt, xt−1, . . .)] (3.17) dla wszystkich s > 0, gdzie MSE oznacza średni błąd prognozy. Literatura poświęcona badaniu przyczynowości i zależności między energią a wzrostem go- spodarczym jest bardzo szeroka. Analizy poświęcone są przypadkom bardzo wielu krajów, np.: Wolde-Rufael62 dla 17 krajów afrykańskich, Squalli63 dla 11 krajów eksportujących ropy naftowej czy Lee i innych64 dla 22 krajów OECD, oraz są zróżnicowane pod względem wykorzystanych metod, np.: przyczynowość Grangera w Kraft i Kraft65, badanie koin- tegracji w Abosedra i Baghestani66, wykorzystanie analizy falkowej w Aslan i innych67, czy model keynesowski w Narayan i innych68. Warto tutaj także przytoczyć artykuł Yu i Choi69, w którym w grupie analizowanych krajów znalazła się też Polska. Obszerna i szczegółowa analiza literatury poświęconej zagadnieniom istnienia oraz kie- runku zależności przyczynowo-skutkowej między energią a PKB przedstawiona jest w Ozturk70. Choć wiedza o występowaniu związku między tymi zmiennymi może być nie- zwykle przydatna dla polityków i decydentów, to nie istnieją jednoznaczne wnioski po- twierdzające uniwersalny jego kierunek. W zależności od użytego zbioru danych, hory- zontu czasowego oraz zastosowanej metody wnioski się różnią. Ozturk71 przytacza bada- Methods, Econometrica, 37(3), s. 424–438. 61Hamilton J.D., 1994, Time Series Analysis, Princeton University Press, Princeton, str. 303. 62Wolde-Rufael Y., 2006, Electricity Consumption and Economic Growth: A Time Series Experience for 17 African Countries, Energy Policy, 34(10), s. 1106-1114. 63Squalli J., 2007, Electricity Consumption and Economic Growth: Bounds and Causality Analyses for OPEC members, Energy Economics, 29(6), s. 1192-1205. 64Lee C.-C., Chang C.-P., Chen P.-F., 2008, Energy Consumption and GDP in Developing Countries: A Cointegrated Panel Analysis, Energy Economics, 27(3), s. 415-427. 65Kraft J., Kraft A., 1978, On the Relationship Between Energy and GNP, Journal of the Energy Development, 3(2), s. 401-403. 66Abosedra S., Baghestani H., 1989, New Evidence on the Causal Relationship between U.S. Energy, Consumption and Gross National Product, Journal of Energy Development, 14, s. 285-292. 67Aslan A., Apergis N., Yildirim S., 2014, Causality between Energy Consumption and GDP in the U.S.: Evidence from Wavelet Analysis, Frontiers in Energy, 8(1), s. 1-8. 68Narayan P.K., Narayan S., Smyth R., 2011, Energy Consumption at Business Cycle Horizons: The Case of the United States, Energy Economics, 33, s. 161-167. 69Yu E.S.H., Choi J.Y., 1985, The Causal Relationship between Energy and GNP: An International Comparison, Journal of Energy and Development, 10, s. 249–272. 70Ozturk I., 2010, A Literature Survey on Energy-Growth Nexus, Energy Policy, 38(1), s. 340-349. 71Ozturk I., 2010, A Literature Survey... 3.5. ZMIENNE OBJAŚNIAJĄCE 65 objaśniającą xt i ma postać yt = µt + k∑ i=1 γi,t + βxt + εt, (3.18) µt = µt−1 + νt−1, (3.19) νt = νt + ζt, (3.20) γi,t = γi,t−1 cos (λi) + γ∗i,t−1 sin (λi) + ωi,t, (3.21) γ∗i,t = −γi,t−1 sin (λi) + γ∗i,t−1 cos (λi) + ω∗i,t, (3.22) gdzie β jest nieznanym współczynnikiem. Estymacja modelu (3.18)-(3.22) jest analogiczna do estymacji modeli przedstawionych w rozdziale 2. Wprowadzenie zmiennej objaśniającej do modelu przekłada się na zmianę w filtrze Kalmana w etapie predykcji, w którym wartość błędu vt|t−1 jest wyznaczana w sposób vt|t−1 = yt −Hαt|t−1 − A′Xt, gdzie Xt jest macierzą zmiennych objaśniających. Rysunek 3.8: Zmienna objaśniająca skonstruowana z komponentu cyklicznego miesięcz- nego PKB Źródło: Opracowanie własne. Wyestymowany w rozdziale 2 model potrójnej sezonowości trygonometrycznej z determi- nistycznym trendem, jak już zostało zauważone, nie bierze pod uwagę wpływu wahań średniookresowych na popyt na energię elektryczną. Uwzględnienie oddziaływania cy- klu koniunkturalnego jest możliwe w takim modelu na przykład przez dodanie do niego zmiennej objaśniającej. Wyestymowane w podrozdziale 3.2 miesięczne PKB zostało użyte jako zmienna objaśniająca w modelach wysokoczęstotliwościowych zapotrzebowania na 66 ROZDZIAŁ 3. WPŁYW CYKLU KONIUNKTURALNEGO Tablica 3.10: Wyniki estymacji modeli zawierających miesięczne PKB jako zmienną ob- jaśniającą Model z podwójną sezonowością Model z potrójną sezonowością + MPKB + MPKB 105σ2 ε 2,692 2,691 1,334 1,333 105σ2 η 3,636 3,636 0 0 105σ2 ω 3,467 3,473 139,138 118,196 β -0,011 -0,012 10−3 logL 2,135126 2,135127 2,225767 2,225769 R2 53,141% 53,142% 79,260% 79,263% Źródło: Opracowanie własne. energię elektryczną. Miesięczne PKB mają częstotliwość miesięczną i jest to najwyższa częstość z jaką, nie budząc kontrowersji, można aproksymować wahania koniunkturalne. W związku z tym, miesięczne wskaźniki zostały dodane do modelu 15-minutowego w taki sposób, że dla każdego miesiąca jego wartość została powtórzona dla wszystkich kwadran- sów. Oznacza to, że wartości x ma postać funkcji schodkowej i jej wartości zmieniają się tylko w momentach zmiany miesiąca. Ilustracją tego jest rysunek 3.8, na którym wyryso- wane zostały znormalizowane i pozbawione średnich szeregi logarytmu zapotrzebowania na energię i zmiennej x dla lat 2011-2012. Z powodu uwzględnienia trendu w modelowanych danych, zmienne objaśniające zostały uwzględnione w wersji stacjonarnej - jako swoje komponenty cykliczne. Tablica 3.10 prezentuje oszacowania modeli z dodatkową zmienną obrazującą koniunkturę dla okresu 2011-2012. Niewielka zmiana wartości współczynnika determinacji oznacza, że dołączenie do modelu zmiennej objaśniającej nie poprawia znacznie dopasowania modelu do danych. Formalne przetestowanie istotności parametru β wykonane zostało za pomocą testu ilo- razu wiarogodności. Wartości statystyk testowych dla modeli z podwójną i potrojną strukturą sezonową wyniosły odpowiednio 0,2 i 0,4, co przy wartości krytycznej równej 3,84 nie daje podstaw81 do odrzucenia hipotezy o nieistotności parametru β. Prawdo- podobną tego przyczyną może być długość próbki, a model ze zmienną aproksymującą wahania średniookresowe należałoby oszacować dla większego zakresu czasowego danych. 81Odpowiadające statystykom wartości p-value były równe odpowiednio 0,6547 i 0,5271. Rozdział 4 Zmienne egzogeniczne W poprzednich rozdziałach wykorzystane zostały jedynie informacje zawarte w szeregu czasowym zapotrzebowania na moce przesyłowe krajowych sieci energetycznych. Oprócz analizy składników cyklu, trendu i sezonowego istotne z punktu widzenia modelowania mogą być informacje zewnętrzne. Najoczywistszym jest uwzględnienie zmiennych pogo- dowych. Warunki atmosferyczne mają wpływ zarówno na wiele decyzji podmiotów. O ile część wahań popytu na energię może być przypisana do stałych efektów kalendarzowych, na przykład do godziny wschodu i zachodu słońca, to inna część, przypisywana warunkom atmosferycznym, dla przykladu zachmurzeniu, nie może mieć deterministycznej postaci. Intuicja ekonomiczna wskazuje na istnienie istotnej zależności między wyborami gospoda- rujących podmiotów a bieżącymi i prognozowanymi warunkami atmosferycznymi. Spędze- nie weekendowego popołudnia w parku czy na działce wydaje się bardziej prawdopodobne podczas pięknej, słonecznej, złotej polskiej jesieni niż w chłodny deszczowy dzień końca września. Dodatkowo, przyjemna temperatura wieczorem może skłaniać do zorganizowa- nia ogniska lub grilla, co także wpływa na schemat zużycia energii. Wydaje się, że sobotnie pochmurne popołudnie spędzone w domu powinno się przełożyć na wzrost poboru energii elektrycznej w stosunku do dnia spędzonego poza domem przy pozostałych warunkach niezmienionych. Wysokie temperatury latem i niskie zimą mogą przekładać się na wzrost zapotrzebowania na wzrost popytu na energię ze względu na wykorzystanie klimatyzacji i ogrzewania. Uzasadnioną wydaje się być w związku z tym hipoteza o wpływie czynników atmosferycznych na popyt na energię elektryczną. W bieżącym rozdziale analizowany jest wpływ zmiennych pogodowych na zapotrzebo- wanie na moce przesyłowe krajowych sieci elektroenergetycznych za pomocą różnych - standardowych i mniej standardowych - modeli. W podrozdziale 4.1 zaprezentowane są dane pogodowe. Podrozdział 4.2 poświęcony jest analizie danych miesięcznych, co spowo- dowane jest potrzebą dalszej weryfikacji hipotezy o wpływie cyklu koniunkturalnego na energię przesyłaną w sieci. W ramach tej analizy uwzględniłem także informację o warun- 67 70 ROZDZIAŁ 4. ZMIENNE EGZOGENICZNE Dane przedstawiają sumaryczny obraz stanu pogody dla wszystkich, dostępnych serwi- sowi ogimet.com stacji meteorologicznych. Zawarte w nich są informacje o temperaturze (maksimum, minimum i mediana), temperaturze punktu rosy, wilgotności, wietrze (kie- runku, intensywności i porywie), ciśnieniu atmosferycznym, wielkości opadów, poziomie zachmurzenia, nasłonecznieniu, widoczności poziomej i pokrywie śnieżnej. Fragment da- nych został przestawiony w tablicy 4.1. Tablica 4.2: Podstawowe statystyki zmiennych pogodowych Zmienna Liczba obs. Średnia Odch. std. Min Max lowclocta 103211 4,34 1,74 0,00 7,70 |∆lowclocta| 103209 0,05 0,09 0,00 3,30 temp_med 103211 8,47 8,31 -20,70 26,50 |∆temp_med| 103209 0,09 0,09 0,00 4,40 visibility 103211 14,90 6,27 2,40 35,00 |∆visibility| 103209 0,14 0,24 0,00 4,30 Źródło: Opracowanie własne. W przedstawionych w dalszej części pracy wykorzystałem medianę temperatury, zachmu- rzenie chmurami piętra niskiego10, widoczności i nasłonecznieniu. W modelach dla danych 15-minutowych za warunki atmosferyczne dla poszczególnych kwadransów użyłem warto- ści zaobserwowanych dla godziny. Jestem świadomy, że jest to pewne przybliżenie, jednak analizy danych pogodowych11 wykazały małą zmienność wartości w horyzoncie 2 godzin. Argumentem jest także to, że dane są wartościami średnimi lub medianami danych uzy- skanych dla 76 stacji meteorologicznych. Dodatkowym uzasadnieniem jest także to, że w polskim klimacie ekstremalne zjawiska pogodowe zdarzają się niezwykle rzadko. 4.2 Analiza miesięczna Konwersja danych wysokoczęstotliwościowych na niskoczęstotliwościowe zazwyczaj polega na zsumowaniu obserwacji z odpowiednich okresów. W przypadku danych o wykorzy- staniu mocy przesyłowych krajowych sieci elektroenergetycznych postępowanie musi być odmienne. Zmienna ta nie jest zmienną przepływu lecz stanu i w związku z tym należy postępować jak na przykład z temperaturą czy poziomem wody w rzece - czyli wybrać jedną z obserwacji z agregowanego okresu. Często wykorzystywana w takich przypadkach „reguła kciuka” polega na wyborze obserwacji z trzeciej środy każdego miesiąca. Postępo- wanie to ma zapewnić lepszą porównywalność danych między miesiącami oraz usunięcie części - niechcianych i zaburzających wnioskowanie - efektów sezonowych. W literaturze 10Na podstawie wstępnych analiz wyciągnąłem wniosek o istotności zachmurzenia chmurami piętra niskiego oraz o nieistotności całkowitego zachmurzenia. 11Odchylenia standardowe pierwszych i drugich różnic szeregów widoczności, zachmurzenia i tempera- tury nie przekraczały wartości 0,3. 4.2. ANALIZA MIESIĘCZNA 71 zabieg ten można najczęściej spotkać spotkać w publikacjach dotyczących testowania hi- potezy o efektywności rynku finansowego. Bigman, Goldfarb i Schechtman12 analizowali efektywność rynku, gdzie między innymi cena spot pochodziła z trzeciej środy miesiąca zakończenia kontraktu. Anderson, Duru i Reeb13 użyli danych z trzeciej środy miesiąca do wyliczenia spreadu akcji w analizie korporacyjnej nieprzejrzystości spółek w Stanach Zjednoczonych. Wybór ten został tam uzasadniony wyeliminowaniem braków danych spo- wodowanych świętami oraz minimalizacją efektu sezonowego dni tygodnia.14Dodatkowego argumentu dostarczają także publikacje Agencji Rynku Energii15, w których prezentowane są dobowe wykresy zapotrzebowania na moce przesyłowe KSE w trzecią środę miesiąca. Jak już zostało wspomniane na początku tej części pracy poświęconej zmiennym egzoge- nicznym, modele przedstawione w tym podrozdziale są modelami dla danych miesięcznych. Postępowanie to ma pozwolić zweryfikować hipotezę o sensowności wykorzystania infor- macji o cyklu koniunkturalnym w modelowaniu popytu na energię elektryczną w szerszym zakresie czasowym, ale przy użyciu klasycznego narzędzia – regresji liniowej. Tablica 4.3 przedstawia oszacowania parametrów modeli regresji liniowej dla danych miesięcznych oraz standardowy i skorygowany współczynnik R2 i R2 a. Wiersze White, Ramsey i BGod- frey zawierają wartości p-value dla testowania hipotez odpowiednio o homoskedastyczności reszt, poprawności formy funkcyjnej oraz braku autokorelacji reszt odpowiedniego rzędu. We wszystkich poniższych przypadkach modelowaną zmienną y jest logarytm kwadran- sowego zapotrzebowania na moce przesyłowe KSE w GWh. Oszacowane zostały modele regresji liniowej postaci: yt = β0 + β1t+ β2miesiac 2 t + . . .+ β12miesiac 12 t + β13coint + β14t ∗ coint+ (4.1) +β15miesiac 1 t ∗ t+ . . .+ β26miesiac 12 t ∗ t+ εt, gdzie εt ∼ NID(0, σ2 ε), a zmienna coin oznacza komponent cykliczny wskaźnika koniunk- tury przedstawionego w rozdziale 3. Przedstawione modele zawierają komponent determi- nistycznej sezonowości modelowanej za pomocą zmiennych zero-jedynkowych oraz deter- ministycznego liniowego trendu (t). Zmienna dla stycznia została pominięta, ze względu na problem współliniowości ze stałą i jako taka jest poziomem bazowym. Modele (1) i (2) z tablicy 4.3 miały pozwolić na wnioskowanie o istotności cyklu koniunkturalnego w modelowaniu zapotrzebowania na moce przesyłowe, jednak ze względu na niskie wartości p-value, w testach na poprawność formy funkcyjnej i homoskedastyczność, wnioskowanie na ich podstawie jest nieuprawnione. Rozszerzenie modelu o interakcję między determini- stycznym trendem liniowym a zmiennymi zero-jedynkowymi dla poszczególnych miesięcy pozwala rozwiązać problemy, występujące wcześniej, z niepoprawną formą funkcyjną i 12Bigman D., Goldfarb D., Schechtman E., 1983, The Journal of Futures Markets, 3(3), s. 321-334. 13Anderson R.C., Duru A., Reeb D.M., 2009, Founders, Heirs, and Corporate Opacity in the United States, Journal of Financial Economics, 92, s. 205–222. 14Tamże, str. 209. 15ARE 2004, Statystyka elektroenergetyki polskiej, Agencja Rynku Energii, Warszawa, str. 44. 72 ROZDZIAŁ 4. ZMIENNE EGZOGENICZNE Tablica 4.3: Miesięczne modele regresji (1) (2) (3) (4) (5) y y y y y t 0,001∗∗∗ 0,001∗∗∗ 0,001∗∗∗ 0,001∗∗∗ 0,001∗∗∗ luty −0,019 −0,019 −0,120 −0,120 −0,118 marzec −0,074∗∗∗ −0,074∗∗∗ −0,321 −0,321 −0,317 kwiecień −0,137∗∗∗ −0,138∗∗∗ −0,327 −0,327 −0,320 maj −0,167∗∗∗ −0,168∗∗∗ −0,586∗ −0,585∗ −0,576∗ czerwiec −0,179∗∗∗ −0,179∗∗∗ −0,698∗∗ −0,697∗∗ −0,686∗∗ lipiec −0,180∗∗∗ −0,181∗∗∗ −0,664∗∗ −0,663∗∗ −0,650∗∗ sierpień −0,231∗∗∗ −0,232∗∗∗ −0,146 −0,145 −0,131 wrzesień −0,153∗∗∗ −0,154∗∗∗ −0,444 −0,443 −0,427 październik −0,083∗∗∗ −0,084∗∗∗ −0,149 −0,147 −0,135 listopad −0,050∗∗ −0,051∗∗ −0,051 −0,049 −0,033 grudzień −0,000 −0,001 0,070 0,072 0,094 coin 0,298∗∗ 0,475 0,297∗∗ t∗coin −0,000 styczeń∗t 0,000 0,000 0,000 luty∗t 0,000 0,000 0,000 marzec∗t 0,000 0,000 0,000 kwiecień∗t 0,000 0,000 0,000 maj∗t 0,001 0,001 0,001 czerwiec∗t 0,001∗ 0,001∗ 0,001∗ lipiec∗t 0,001∗ 0,001∗ 0,001∗ sierpień∗t −0,000 −0,000 −0,000 wrzesień∗t 0,001 0,001 0,000 październik∗t 0,000 0,000 0,000 listopad∗t 0,000 −0,000 −0,000 grudzień∗t −0,000 −0,000 −0,000 Stała 9,161∗∗∗ 9,196∗∗∗ 9,389∗∗∗ 9,386∗∗∗ 9,341∗∗∗ R2 0,845 0,855 0,875 0,875 0,865 R2 a 0,830 0,840 0,848 0,849 0,838 N 141,000 141,000 141,000 141,000 141,000 Ramsey 0,003 0,037 0,360 0,304 0,011 White 0,002 0,025 0,906 0,938 0,883 BGodfrey L1 0,022 0,114 0,207 0,206 0,038 L2 0,007 0,063 0,099 0,986 0,010 L12 0,323 0,622 0,155 0,153 0,050 ∗ p < 0,05, ∗∗ p < 0,01, ∗∗∗ p < 0,001 Źródło: Opracowanie własne. 4.3. ANALIZA DANYCH DZIENNYCH 75 4.3 Analiza danych dziennych Rozważania z poprzedniego rozdziału dają przesłanki do uwzględnienia informacji o cyklu koniunkturalnym i temperaturze w modelowaniu zapotrzebowania na moce przesyłowe sieci energetycznych. Zwiększenie częstotliwości, z miesięcznej na dzienną, ma na celu łatwiejsze wychwycenie występujących w danych zależności oraz w mniejszym stopniu argumentowane jest szybkością estymacji modeli z mniejszą liczbą obserwacji. Podstawą analizy było 4290 obserwacji z godziny 11:15 dla wszystkich dni kalendarzowych ze zbioru danych. Tablica 4.5 zawiera oszacowania trzech modeli regresji dla danych dziennych. Pierwszy z nich jest próbą jak najprostszego modelowania tych obserwacji, ale niestety hipoteza o poprawności formy funkcyjnej jest silnie odrzucana (p-value testu Ramseya równe 0). Intuicyjnie fakt ten można zinterpretować w ten sposób, że wzrost częstotliwości danych powinien przekładać na wzrost skomplikowania modelu. Model (2) jest odpowiednikiem modelu (2) z tablicy 4.4 estymowanym dla dziennych danych, jednak z analogicznego powodu co model (1) jest on odrzucany. Dołączenie interakcji między dniem a miesiącem, ze względu na odmienny wpływ dni w poszczególnych miesiącach, pozwoliło uzyskać model (3) o poprawnej formie funkcyjnej17. Warte podkreślenia jest potwierdzenie się przypuszczenia z podrozdziału 4.2 o tym, że pozatemperaturowe zmienne pogodowe mogą okazać się istotne dla danych o wyższej częstotliwości. Zrozumiałe są znaki oszacowań przy poziomie zachmurzenia (lowclocta) i nasłonecznieniem (sun). Intuicyjne okazało się ponadto oszacowanie przy widoczności (visibility), które należy interpretować jako spadek zapotrzebowania na moce przesyłowe KSE średnio o 0,1% wraz z jednostkowym wzrostem poziomu widoczności. Zwraca uwagę, zarówno w tablicy 4.4 jak i 4.5, znak oszacowania parametru przy medianie temperatury. W modelu (4) wartość oszacowania -0,005 oznacza, że wzrost temperatury o jeden stopień Celsjusza przełoży się średnio na spadek zapotrzebowania na moce przesy- łowe o 0,5%. Wniosek wydaje się być uzasadniony podczas jesieni, zimy i wiosny, jednak latem sytuacja może być odmienna. Gdy tylko temperatura przewyższa pewien próg, zaobserwować można wzrost poboru energii elektrycznej związanej najprawdopodobniej z działaniem urządzeń klimatyzacyjnych. W związku z tym, uzasadnione wydaje się być odwołanie do bardziej zaawansowanych technik modelowania nieliniowego wpływu tem- peratury na popyt na energię elektryczną. 17W modelu tym spełnione jest też założenie o braku autokorelacji reszt. 76 ROZDZIAŁ 4. ZMIENNE EGZOGENICZNE Tablica 4.5: Modele regresji dla danych dziennych (1) (2) (3) y y y t 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ coin 0,484∗∗∗ 0,494∗∗∗ 0,494∗∗∗ temp_med −0,004∗∗∗ −0,004∗∗∗ −0,005∗∗∗ lowclocta 0,004∗∗∗ 0,003∗∗∗ 0,003∗∗∗ sun 0,002∗∗∗ 0,002∗∗∗ 0,002∗∗∗ visibility −0,000 −0,001∗∗∗ −0,001∗∗∗ poniedziałek 0,236∗∗∗ 0,236∗∗∗ 0,248∗∗∗ wtorek 0,254∗∗∗ 0,254∗∗∗ 0,244∗∗∗ środa 0,260∗∗∗ 0,260∗∗∗ 0,257∗∗∗ czwartek 0,251∗∗∗ 0,251∗∗∗ 0,244∗∗∗ piątek 0,250∗∗∗ 0,250∗∗∗ 0,239∗∗∗ sobota 0,157∗∗∗ 0,157∗∗∗ 0,144∗∗∗ luty 0,002 −0,014 −0,015 marzec −0,032∗∗∗ −0,049∗∗∗ −0,055∗∗∗ kwiecień −0,084∗∗∗ −0,110∗∗∗ −0,125∗∗∗ maj −0,111∗∗∗ −0,160∗∗∗ −0,157∗∗∗ czerwiec −0,087∗∗∗ −0,150∗∗∗ −0,157∗∗∗ lipiec −0,072∗∗∗ −0,136∗∗∗ −0,152∗∗∗ sierpień −0,081∗∗∗ −0,139∗∗∗ −0,148∗∗∗ wrzesień −0,065∗∗∗ −0,111∗∗∗ −0,120∗∗∗ październik −0,026∗∗∗ −0,048∗∗∗ −0,057∗∗∗ listopad −0,023∗∗∗ −0,035∗∗∗ −0,033∗ grudzień −0,021∗∗∗ −0,015 0,004 luty∗t 0,000∗ 0,000∗ marzec∗t 0,000∗ 0,000∗ kwiecień∗t 0,000∗∗∗ 0,000∗∗∗ maj∗t 0,000∗∗∗ 0,000∗∗∗ czerwiec∗t 0,000∗∗∗ 0,000∗∗∗ lipiec∗t 0,000∗∗∗ 0,000∗∗∗ sierpień∗t 0,000∗∗∗ 0,000∗∗∗ wrzesień∗t 0,000∗∗∗ 0,000∗∗∗ październik∗t 0,000∗∗ 0,000∗∗ listopad∗t 0,000 0,000 grudzień∗t −0,000 −0,000 poniedziałek & luty −0,007 ... ... sobota & grudzień −0,011 Stała 9,586∗∗∗ 9,627∗∗∗ 9,631∗∗∗ R2 0,805 0,815 0,821 R2 a 0,804 0,814 0,816 N 4290,000 4290,000 4290,000 Ramsey 0,000 0,003 0,050 ∗ p < 0,05, ∗∗ p < 0,01, ∗∗∗ p < 0,001 Źródło: Opracowanie własne. 4.4. NIELINIOWE ZALEŻNOŚCI 77 4.4 Nieliniowe zależności Rozgraniczenie w literaturze przedmiotu na stopniodni grzania (heating degree day) i stop- niodni chłodzenia (cooling degree day) świadczy o nieliniowym przełożeniu tempretury na zapotrzebowanie na energię. Nieliniowa zależność między tymi zmiennymi w literaturze pojawia się najczęściej w postaci V-kształtnej. Engle i inni18, Harvey i Koopman19, Valor, Meneu i Caselles20 oraz Dordonnat i inni21 uzyskiwali dla danych ze Stanów Zjednoczo- nych, Hiszpanii i Francji taką właśnie zależność. Na tej podstawie można wysunąć hipotezę o tym, że także w przypadku danych polskich taką V-kształtnosć zależności powinno się zaobserwować. Wzrost zapotrzebowania na energię elektryczną w okresie letnim związany jest z wysokimi temperaturami, które po- wodują wzrost wykorzystania urządzeń chłodzących, w tym klimatyzacji, przez co wzrasta zapotrzebowanie na energię. Argumentem przeciwnym jest niska liczba domów i miesz- kań wyposażonych w urządzenia klimatyzacyjne22. Według statystyki Głównego Urzędu Statystycznego urządzenia klimatyzacyjne znajdowały się w 0,35% gospodarstw domo- wych23, podczas gdy w dane amerykańskie wskazują na analogiczny odsetek wynoszący prawie 83%24. Podejrzewać można jednak, że w biurach, urzędach i centrach handlowych w Polsce wykorzystuje się dużą liczbę urządzeń klimatyzacyjnych, więc w dalszej części rozważam hipotezę o nieliniowości relacji między temperaturą a zapotrzebowaniem na energię. Nieliniowa zależność może być analizowana w modelu semiparametrycznym lub poprzez liniowe przybliżenie w modelu regresji. Zaprezentuję najpierw wyniki dwóch modeli semi- parametrycznych, po czym zaproponuję sposób modelowania zależności między popytem na energię elektryczną a temperaturą w liniowym modelu szacowanym metodą najmniej- szych kwadratów. Model regresji semiparametrycznej ma postać yi = β0 + β1x1,i + . . .+ βKxK,i + f(zi) + εi, (4.2) dla i = 1, 2, . . . , N , gdzie y jest zmienną objaśnianą, x1, . . . , xK są zmiennymi niezależ- nymi, ε jest zaburzeniem z zerową średnią i stałą wariancją σ2 ε , natomiast z jest zmienną 18Engle R.F., Granger C.W.J., Rice J., Weiss A., 1986, Semiparametric Estimates... 19Harvey A.C., Koopman S.J., 1993, Forecasting Hourly Electricity Demand... 20Valor E., Meneu V., Caselles V., 2001, Daily Air Temperature... 21Dordonnat V., Koopman S.J., Ooms M., Dessertaine A., 2008, An Hourly Periodic State Space Model... 22Przesłanką potwierdzającą hipotezę mógłaby być statystyka urządzeń klimatyzacyjnych w biurach, urzędach i lokalach handlowo-usługowych. Dane dla Stanów Zjednoczonych publikowane są przez Energy Information Administration(U.S. Energy Information Administration, 2009, Commercial Building Energy Consumption Survey), natomiast ich odpowiednik dla Polski nie istnieje. 23GUS, 2014, Zużycie energii w gospodarstwach domowych w 2012 r., Główny Urząd Statystyczny, Warszawa, str. 36. 24U.S. Energy Information Administration, 2009, Residential Energy Consumption Survey, tablica HC 7.1 Air Conditioning in U.S. Homes, by Housing Unit Type: Na 113,6 mln gospodarstw domowych dane wskazują, że 94 mln posiadało sprawne i używało urządzenia klimatyzacji. 80 ROZDZIAŁ 4. ZMIENNE EGZOGENICZNE korzystaniem procedur odpowiednio Robinsona i Yatchewa, natomiast (3) i (4) modelo- wane w sposób parametryczny przy użyciu metody najmniejszych kwadratów. Model (3) jest modelem regresji z funkcją łamaną, a (4) można obrazowo ująć jako „kij hokejowy”34 lub małe „l” „pisane”. Niech pt będzie sumą wszystkich zmiennych pozatemperaturowych wpływających na zapotrzebowanie na moce przesyłowe pt = β0 + β1t+ β2coint + β3dzient + β4miesiact + β5dzient ∗miesiact+ + β6t ∗miesiact + β7lowcloctat + β8sunt + β9visibilityt, (4.13) wówczas model (3) można zapisać jako yt = pt+β10temp_medt+β111(temp_medt>16)+β121(temp_medt>16)∗temp_medt+εt, (4.14) a model (4) w postaci yt = pt + β10temp_medt + 1(14<temp_medt<18) ( β11 + β12(temp_medt − 16)2 ) + + 1(temp_medt>18) (β13 + β14temp_medt) + εt. (4.15) Rysunek 4.1: Nieliniowy wpływ temperatury dla danych dziennych 9. 4 9. 6 9. 8 10 10 .2 −20 −10 0 10 20 30 temp_med Źródło: Opracowanie własne. 34Weron R., 2006, Modeling and Forecasting..., str. 74 określił w ten sposób zależność między tempa- raturą a zapotrzebowaniem na energię elektryczną. 4.4. NIELINIOWE ZALEŻNOŚCI 81 Tablica 4.6: Dzienne modele z nieliniową zależnością między energią a temperaturą (1) (2) (3) (4) y y y y t 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ coin 0,516∗∗∗ 0,433∗∗∗ 0,509∗∗∗ 0,511∗∗∗ lowclocta 0,006∗∗∗ 0,005∗∗∗ 0,005∗∗∗ 0,005∗∗∗ sun 0,001∗ 0,001 0,001∗ 0,001∗ visibility −0,000 −0,001∗∗∗ −0,000 −0,000 poniedziałek 0,245∗∗∗ 0,251∗∗∗ 0,248∗∗∗ 0,249∗∗∗ wtorek 0,239∗∗∗ 0,238∗∗∗ 0,246∗∗∗ 0,246∗∗∗ środa 0,256∗∗∗ 0,266∗∗∗ 0,260∗∗∗ 0,260∗∗∗ czwartek 0,241∗∗∗ 0,234∗∗∗ 0,246∗∗∗ 0,246∗∗∗ piątek 0,235∗∗∗ 0,236∗∗∗ 0,239∗∗∗ 0,239∗∗∗ sobota 0,142∗∗∗ 0,141∗∗∗ 0,144∗∗∗ 0,144∗∗∗ luty −0,012 −0,022 −0,011 −0,010 marzec −0,047∗∗∗ −0,048∗∗ −0,044∗∗ −0,043∗∗ ... ... ... ... ... grudzień 0,004 0,004 0,008 0,009 poniedziałek & luty −0,006 0,007 −0,007 −0,008 poniedziałek & marzec −0,010 −0,010 −0,014 −0,014 poniedziałek & kwiecień −0,052∗∗ −0,058∗∗ −0,055∗∗∗ −0,056∗∗∗ ... ... ... ... ... sobota & grudzień −0,007 0,016 −0,010 −0,010 luty*t 0,000 0,000 0,000 0,000 marzec*t 0,000∗∗ 0,000∗ 0,000∗ 0,000∗ kwiecień*t 0,000∗∗∗ 0,000∗∗ 0,000∗∗∗ 0,000∗∗∗ maj*t 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ czerwiec*t 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ lipiec*t 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ sierpień*t 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ wrzesień*t 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ październik*t 0,000∗∗∗ 0,000∗ 0,000∗∗∗ 0,000∗∗∗ listopad*t 0,000∗ 0,000 0,000∗ 0,000∗ grudzień*t −0,000 −0,000 −0,000 −0,000 temp_med −0,006∗∗∗ 1(temp_medt>16) −0,189∗∗∗ 1(temp_medt>16)temp_med 0,012∗∗∗ 1(temp_medt<14)temp_med −0,006∗∗∗ 1(14≤temp_medt≤18) −0,088∗∗∗ 1(14≤temp_medt≤18)(temp_med-16)2 0,002 1(temp_medt>18) −0,213∗∗∗ 1(temp_medt>18)temp_med 0,007∗∗∗ Stała 9,614∗∗∗ 9,613∗∗∗ R2 0,785 0,752 0,831 0,831 R2 a 0,780 0,746 0,827 0,827 N 4290,000 4289,000 4290,000 4290,000 ∗ p < 0,05, ∗∗ p < 0,01, ∗∗∗ p < 0,001 Źródło: Opracowanie własne. 82 ROZDZIAŁ 4. ZMIENNE EGZOGENICZNE Zaprezentowane w tablicy 4.6 modele, będące rozszerzeniem (4) z tablicy 4.5, potwier- dzają hipotezę o istotności cyklu koniunkturalnego oraz nieliniowym wpływie tempera- tury w przypadku danych dziennych35. Nieliniowość ta jest zaprezentowana na rysunku 4.1, który prezentuje wynik estymacji modelu (4.2) metodą Robinsona. Pokazany zo- stał tylko wykres uzyskany z metody Robinsona ze względu na niemalże całkowite jego podobieństwo do relacji wyestymowanej metodą Yatchewa. Oszacowania modelu (3) i (4) potwierdzają hipotezę o nieliniowości zależności między temperaturą a popytem na energię, a ich wizualizacje zamieszczone są na rysunku 4.2. Potwierdzone zostało ponadto przypuszczenie o V-kształtności relacji, ponieważ dla chłod- nych dni spadek temperatury o 1 stopień Celsjusza powoduje średnio spadek zapotrzebo- wania na energię o 0,6%, natomiast w ciepłe dni oczekujemy, że wzrostowi temperatury o 1◦C towarzyszyć będzie wzrost zapotrzebowania na energię o 0,6%. Bardzo podobne wnioski mogą być wyciągnięte na podstawie oszacowań modelu (4), w którym charakter zależności zmienia się bardziej gładko, to jest w otoczeniu temperatury 16◦C zgodnie z funkcją kwadratową. Obserwowany jest analogiczny wzrost podczas dni grzewczych, pod- czas gdy dla temperatur przekraczających 18◦C wzrost temperatury prowadzi do wyższego wzrostu zapotrzebowania na moce przesyłowe, to jest o 0,7% przy 0,6% wcześniej. Za- uważyć należy, że na podstawie oszacowań uzyskanych dla próby danych dziennych nie można wnioskować, który ze sposobów parametrycznego uwzględniania temperatury jest lepszy, dlatego, że obydwa modele mają równe wartości współczynników determinacji. Rysunek 4.2: Wykres oszacowanych parametrycznych nielinowych zależności −20 −10 0 10 20 30 9.5 9.55 9.6 9.65 9.7 9.75 krzywa lamana kij hokejowy Źródło: Opracowanie własne. 35Nie jest to widoczne bezpośrednio w tablicy 4.6 na podstawie oszacowań, ponieważ relacja była estymowana nieparametrycznie. 4.5. WPŁYW DNI WOLNYCH OD PRACY 85 Zastosowanie zmiennych dla pierwszego i szóstego stycznia, pierwszego i trzeciego maja, jedenastego listopada i długich weekendów bez wątpienia jest kluczowe w modelowaniu dziennego zapotrzebowania na moce przesyłowe KSE. Wykorzystanie pozostałych infor- macji o dniach wolnych ustawowo zostanie sprawdzone także w przypadku danych o wyż- szej częstotliwości. Wyniki uzyskane z mdoeli parametrycznych i semiparametrycznych wskazują na nieli- niowy charakter zależności między temperaturą a zapotrzebowaniem na energię elek- tryczną. W związku z tym pierwsza główna hipoteza badawcza mojej pracy musi zostać zweryfikwoana pozytywnie. Na podstawie oszacowań modeli z tablicy 4.8 stwierdzam, że pod względem skorygowanego R2 wybrać należałoby model (2) zawierający zmienne zero-jedynkowe dla wszystkich dni ustawowo wolnych od pracy. Przy przyjęciu za kryterium wyboru modelu spełnienia przez model wymogu poprawności formy funkcyjnej wybrać należałoby model (4) lub (5). Wybór najlepszego modelu do prognozowania wartości zapotrzebowania na moce przesyłowe krajowego systemu elektroenergetycznego spośród (3)-(5) nie jest oczywisty. Z jednej strony modele (4) i (5) spełniają założenie o poprawności formy funkcyjnej, z drugiej strony model (2) charakteryzuje się wyższym współczynnikiem R2 dzięki czemu może pozwalać wyznaczać dokładniejsze prognozy. 86 ROZDZIAŁ 4. ZMIENNE EGZOGENICZNE Ta bl ic a 4. 7: D ni w ol ne od pr ac y id łu gi e w ee ke nd y w la ta ch 20 02 -2 01 4 N ow y Tr ze ch Zi el on e B oż e W ni eb ow z. W sz ys t. Św ię to B oż e D łu gi e R ok R ok K ró li W ie lk an oc M aj ów ka Św ią tk i C ia ło N M P Św ię t. N ie po dl . N ar od ze ni e W ee ke nd y 20 02 1 I - 31 II I 1 IV 1 3 V 19 V 30 V 15 V II I 1 X I 11 X I 25 -2 6 X II 2V 31 V 16 V II I 20 03 1 I - 20 -2 1 IV 1 3 V 8 V I 19 V I 15 V II I 1 X I 11 X I 25 -2 6 X II 2V 20 V I 10 X I 20 04 1 I - 11 -1 2 IV 1 3 V 30 V 19 V I 15 V II I 1 X I 11 X I 25 -2 6 X II 2I I 12 X I 20 05 1 I - 27 -2 8 II I 1 3 V 15 V 26 V 15 V II I 1 X I 11 X I 25 -2 6 X II 2V 27 V 31 X 20 06 1 I - 16 -1 7 IV 1 3 V 4 V 15 V 15 V II I 1 X I 11 X I 25 -2 6 X II 2V 5V 14 V II I 20 07 1 I - 8- 9 IV 1 3 V 27 V 7 V I 15 V II I 1 X I 11 X I 25 -2 6 X II 30 IV -4 V 8V I 2X I 31 X II 20 08 1 I - 23 -2 4 II I 1 3 V 11 V 22 V 15 V II I 1 X I 11 X I 25 -2 6 X II 2V 23 V 10 X I 20 09 1 I - 12 -1 3 IV 1 3 V 31 V 11 V I 15 V II I 1 X I 11 X I 25 -2 6 X II 2I 12 V I 20 10 1 I - 4- 5 IV 1 3 V 23 V 3 V I 15 V II I 1 X I 11 X I 25 -2 6 X II 12 X I 20 11 1 I 6 I 24 -2 5 IV 1 3 V 12 V 23 V 15 V II I 1 X I 11 X I 25 -2 6 X II 7I 2V 13 V 31 X 20 12 1 I 6 I 8- 9 IV 1 3 V 27 V 7 V I 15 V II I 1 X I 11 X I 25 -2 6 X II 30 IV -4 V 8V I 2X I 31 X II 20 13 1 I 6 I 31 II I 1 IV 1 3 V 19 V 30 V 15 V II I 1 X I 11 X I 25 -2 6 X II 29 IV -2 V 31 V 16 V II I 27 X II 20 14 1 I 6 I 20 -2 1 IV 1 3 V 8 V I 19 V I 15 V II I 1 X I 11 X I 25 -2 6 X II 2V 20 V I 10 X I Źr ód ło : O pr ac ow an ie w ła sn e. 4.5. WPŁYW DNI WOLNYCH OD PRACY 87 Tablica 4.8: Dzienne modele z nieliniową zależnością między energią a temperaturą i czasem wolnym (1) (2) (3) (4) (5) y y y y y t 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ 0,000∗∗∗ coin 0,509∗∗∗ 0,516∗∗∗ 0,512∗∗∗ 0,508∗∗∗ 0,510∗∗∗ temp_med −0,006∗∗∗ −0,006∗∗∗ −0,006∗∗∗ −0,006∗∗∗ −0,006∗∗∗ lowclocta 0,005∗∗∗ 0,004∗∗∗ 0,004∗∗∗ 0,005∗∗∗ 0,005∗∗∗ sun 0,001∗ 0,001∗∗∗ 0,001∗ 0,001∗ 0,001∗∗ visibility −0,000 −0,000 −0,000 −0,000 −0,000 poniedziałek 0,248∗∗∗ 0,257∗∗∗ 0,256∗∗∗ 0,257∗∗∗ 0,257∗∗∗ wtorek 0,246∗∗∗ 0,263∗∗∗ 0,262∗∗∗ 0,262∗∗∗ 0,263∗∗∗ środa 0,260∗∗∗ 0,268∗∗∗ 0,268∗∗∗ 0,268∗∗∗ 0,268∗∗∗ czwartek 0,246∗∗∗ 0,268∗∗∗ 0,267∗∗∗ 0,267∗∗∗ 0,267∗∗∗ piątek 0,239∗∗∗ 0,260∗∗∗ 0,259∗∗∗ 0,259∗∗∗ 0,259∗∗∗ sobota 0,144∗∗∗ 0,144∗∗∗ 0,144∗∗∗ 0,144∗∗∗ 0,144∗∗∗ luty −0,011 −0,008 −0,008 −0,007 −0,007 marzec −0,044∗∗ −0,017∗ −0,042∗∗∗ −0,041∗∗∗ −0,041∗∗∗ kwiecień −0,103∗∗∗ −0,057∗∗∗ −0,104∗∗∗ −0,101∗∗∗ −0,102∗∗∗ maj −0,129∗∗∗ −0,127∗∗∗ −0,130∗∗∗ −0,127∗∗∗ −0,127∗∗∗ czerwiec −0,136∗∗∗ −0,134∗∗∗ −0,136∗∗∗ −0,134∗∗∗ −0,133∗∗∗ lipiec −0,154∗∗∗ −0,149∗∗∗ −0,152∗∗∗ −0,150∗∗∗ −0,150∗∗∗ sierpień −0,141∗∗∗ −0,136∗∗∗ −0,140∗∗∗ −0,138∗∗∗ −0,137∗∗∗ wrzesień −0,094∗∗∗ −0,092∗∗∗ −0,095∗∗∗ −0,092∗∗∗ −0,092∗∗∗ październik −0,037∗∗ −0,034∗∗∗ −0,036∗∗∗ −0,034∗∗ −0,034∗∗ listopad −0,022 −0,017∗ −0,018 −0,017 −0,017 grudzień 0,008 0,013 0,010 0,011 0,011 ... ... ... ... ... ... grudzień∗t −0,000 −0,000∗ −0,000 −0,000 −0,000 poniedziałek & luty −0,007 −0,015 −0,015 −0,016 −0,016 poniedziałek & marzec −0,014 −0,033∗∗∗ −0,022 −0,022 −0,022 poniedziałek & kwiecień −0,055∗∗∗ −0,047∗∗∗ −0,058∗∗∗ −0,058∗∗∗ −0,058∗∗∗ ... ... ... ... ... ... sobota & grudzień −0,010 −0,010 −0,010 −0,010 −0,010 1(temp_medt>16) −0,189∗∗∗ −0,169∗∗∗ −0,166∗∗∗ −0,171∗∗∗ −0,180∗∗∗ 1(temp_medt>16)temp_med 0,012∗∗∗ 0,011∗∗∗ 0,011∗∗∗ 0,011∗∗∗ 0,012∗∗∗ styczeń 1 −0,444∗∗∗ −0,444∗∗∗ −0,443∗∗∗ −0,443∗∗∗ styczeń 6 −0,268∗∗∗ −0,268∗∗∗ −0,267∗∗∗ −0,267∗∗∗ maj 1 −0,229∗∗∗ −0,229∗∗∗ −0,229∗∗∗ −0,222∗∗∗ maj 3 −0,290∗∗∗ −0,289∗∗∗ −0,289∗∗∗ −0,282∗∗∗ sierpień 15 −0,281∗∗∗ −0,281∗∗∗ listopad 1 −0,401∗∗∗ −0,401∗∗∗ listopad 11 −0,234∗∗∗ −0,234∗∗∗ −0,220∗∗∗ −0,220∗∗∗ grudzień 25 −0,392∗∗∗ grudzień 26 −0,381∗∗∗ Wielkanoc −0,294∗∗∗ Zielone Świątki −0,000 Boże Ciało −0,243∗∗∗ −0,244∗∗∗ −0,244∗∗∗ Długi weekend −0,114∗∗∗ −0,108∗∗∗ −0,106∗∗∗ −0,104∗∗∗ Stała 9,614∗∗∗ 9,615∗∗∗ 9,618∗∗∗ 9,615∗∗∗ 9,613∗∗∗ R2 0,831 0,944 0,899 0,876 0,869 R2 a 0,827 0,943 0,896 0,873 0,865 N 4290,000 4290,000 4290,000 4290,000 4290,000 Ramsey 0,021 0,000 0,012 0,031 0,033 ∗ p < 0,05, ∗∗ p < 0,01, ∗∗∗ p < 0,001 Źródło: Opracowanie własne. 90 ROZDZIAŁ 5. PROGNOZOWANIE POPYTU NA ENERGIĘ ELEKTRYCZNĄ żenia estymatorów i ich ewentualny wpływ na prognozy. Z teorii ekonometrii wynika, że spełnienie założenia o poprawności formy funkcyjnej jest jednoznaczne z nieobciążonością i zgodnością estymatorów, natomiast pozytywne zweryfikowanie hipotez o braku autokore- lacji i heteroskedastyczności pozwala na wnioskowanie statystyczne. Poprawność modelu, rozumiana jako pozytywne zweryfikowanie wszystkich wymaganych postulatów, uzasad- nia prognozy uzyskiwane z modelu, jednak w określonych przypadkach preferowane może być użycie modelu niespełniającego wszystkich założeń teoretycznych. Pojawia się tutaj kolejny problem wyboru między modelem z pozytywnie zweryfikowanymi własnościami a modelem umożliwiającym szybsze lub ekonomiczniejsze prognozy. Powstaje czysto teore- tyczne pytanie o to, czy zawsze model, z którego można uzyskać lepsze prognozy, będzie modelem spełniającym wszystkie założenia. Czy przy prognozowaniu popytu na energię elektryczną można rozluźnić wymóg pozytywnej diagnostyki modelu, jeśli efektem tego byłyby dokładniejsze prognozy? Choć pierwszym wyborem ekonometryka powinna być poprawność statystyczna i ekonometryczna, to podać można przykłady praktycznych roz- wiązań kompromisowych. Dla przykładu, w statystyce małych obszarów2 dopuszcenie obciążoności estymatorów przekłada się na uzyskiwanie oszacowań ze znacząco niższą wa- riancją. W związku z powyższym w tym rozdziale wykorzystane zostaną także metody naiwne, które będą punktem odniesienia do przedstawionych modeli. Ostatnim, choć nie najmniej ważnym, wątkiem jest kwestia modelu używanego do modelo- wania i prognozowania. Pierwszym wyborem zawsze jest model regresji liniowej szacowany metodą najmniejszych kwadratów. Sposób ten umożliwia modelowanie deterministycz- nych zależności, co nie powinno być zawsze postrzegane jako wada. Umiejętność opisania danego zjawiska za pomocą komponentów deterministycznych posiada tę zaletę, że więk- szą część zmienności będziemy w stanie przewidzieć. Wydaje się, że optymalnym roz- wiązaniem byłoby posiadanie modelu regresji liniowej spełniającego wszystkie założenia i posiadającego wysoki współczynnik determinacji liniowej. Niestety, obserwowane zjawiska ekonomiczne charakteryzują się nie deterministycznym, a stochastycznym charakterem. Analiza komponentu sezonowego 15-minutowego szeregu zapotrzebowania na energię elek- tryczną w Polsce, przedstawiona w rozdziale 2, potwierdza stochastyczny charakter tego komponentu, natomiast na podstawie rozdziału pierwszego stwierdzić należy obecność stochastycznego trendu. Wydaje się, że argumenty o stochastycznym charakterze trendu i sezonowości wystarczają do stwierdzenia niemożliwości znalezienia poprawnego modelu regresji liniowej dla 15-minutowych danych Dwie tury zawodów prognostycznych przeprowadzonych w latach 1990-1992 opisane przez 2Statystyka małych obszarów (small area estimation) jest dziedziną statystyki zajmującą się szaco- waniem parametrów dla małych podpopulacji. Mały obszar jest definiowany jako część populacji, która jest w niskim stopniu reprezentowana w próbie. Używanie standardowych narzędzi statystycznych do wyznaczenia charakterystyk małego obszaru wiąże się z dużą wariancją oszacowań. Szczegóły na temat statystyki małych obszarów znaleźć można w Rao N.K.J., 2003, Small Area Estimation, John Wiley & Sons, s. 344. W Polsce głównym ośrodkiem zajmującym się statystyką małych obszarów jest Urząd Statystyczny w Poznaniu i Katedra Statystyki Akademii Ekonomicznej w Poznaniu. 5.1. PROJEKT EKSPERYMENTU 91 Ramanathana i innych3 stanowiły dla mnie podstawę wyboru modeli do prognozowania. Gospodarzem zawodów była spółka energetyczna Puget Sound Power and Light. Udział w zawodach wzięło 11 zespołów badaczy praktyków i teoretyków prognozowania zapotrze- bowania na energię elektryczną. Prognozy były wyliczane o 8:00 dla wszystkich godzin na- stępnego dnia roboczego, natomiast w piątek prognozowano godzinowe zapotrzebowanie na energię elektryczną dla 3 kolejnych dni. Jak zauważyli Ramanathan i inni4 specyfika- cje niektórych modeli nie były upublicznione, natomiast na podstawie wyników zawodów wynika, że zaproponowany przez nich model (EGRV) był bardzo konkurencyjny dla pro- gnoz dla dni roboczych. Model zmiennych w czasie okresowych splajnów kubicznych zaproponowany przez Harvey’a i Koopmana5 w zauważalnej liczbie przypadków okazał się najlepszy. Dla dni weekendowych bardzo dobrze spisywał się model spółki Puget. W aneksie zamieściłem tablicę ze szczegółowymi wynikami zawodów prognostycznych. Ze względu na brak informacji o części modeli w wyborze najlepszego mdoelu dla polskich dancyh wykorzystałem różne wersje modeli EGRV, ARIMA, modelu okresowych splajnów kubicznych. Proponowanie modelu prognostycznego opiera się także na porównaniu jego rezultatów z prognozami naiwnymi. W dalszej części tego rozdziału przedstawię szczegóły eksperymentu, na podstawie którego dokonałem wyboru najlepszego modelu. Podrozdziały 5.1-5.5 zawierają specyfikacje mo- deli, a prognozy z nich uzyskane porównałem w podrozdziale 5.6. Ostatni podrozdział 5.7 zawiera eksperyment za pomocą którego zweryfikowałem drugą główną hipotezę badawczą o przydatności informacji o cyklu koniunkturalnym do prognozowania średniookresowego. 5.1 Projekt eksperymentu Przeprowadzenie eksperymentu porównującego wymaga ustalenia horyzontu prognoz. Li- teratura przedmiotu, na przykład Weron6, dzieli prognozy zapotrzebowania na energię elektryczną na krótko-, średnio- i długo-okresowe. Weron7 zauważył także, że przypisy- wanie konkretnej długości horyzontu prognozy w literaturze nie jest jednoznaczne. Zde- cydowałem, że w eksperymencie zostaną wykorzystane prognozy na 1, 2 i 3 dni jako prognozy krótkookresowe, oraz prognozy na 3 miesiące służące do porównania własności przewidywań modeli w średnim okresie. Analizy i prognozy długookresowe nie zostały przeanalizowane, ponieważ środek ciężkości analizy musiałby się przesunąć z danych o wysokiej częstotliwości w stronę danych o niższych częstościach, prawdopodobnie kwar- 3Ramanathan R., Engle R., Granger C.J.V., Vahid-Araghi F., Brace C., 1997, Short-run Forecast of Electricity Loads and Peaks, International Journal of Forecasting, 13, s. 161-174. 4Ramanathan R., Engle R., Granger C.J.V., Vahid-Araghi F., Brace C., 1997, Short-run Forecast..., str. 168. 5Harvey A.C., Koopman S.J., 1993, Forecasting Hourly Electricity Demand Using Time-Varying Spli- nes, Journal of the American Statitical Association, 88(424), s. 1228-1236. 6Weron R., Modelling and ..., str. 67 7Weron R., Modelling and ..., str. 67 92 ROZDZIAŁ 5. PROGNOZOWANIE POPYTU NA ENERGIĘ ELEKTRYCZNĄ talnych i rocznych. Wydaje się, że najlepszym kryterium do porównania modeli prognostycznych jest spraw- dzenie trafności uzyskiwanych z nich prognoz. Przeprowadziłem w związku z tym ekspe- ryment mający na celu porównanie zaobserwowanych wartości zapotrzebowania na moce przesyłowe krajowych sieci elektroenergetycznych z prognozami porównywanych modeli. Porównanie prognoz krótkookresowych polegało na oszacowaniu modelu na próbie od stycznia 2002 do 10 października 2013 w pierwszym kroku i wyznaczeniu prognoz na trzy kolejne dni. W każdym kolejnym kroku okres, na którym szacowany był model, zwięk- szany był o kolejny dzień. W ten sposób dla każdego modelu przeprowadzonych zostało 445 iteracji, aż do 30 grudnia 2014. Jak już zostało wspomniane, do porównania modeli wybrane zostały: metoda dnia po- dobnego, modele regresji, modele ARIMA i model okresowych splajnów kubicznych. Ex ante nie wiadome było, która możliwa specyfikacja każdego modelu będzie najlepsza. Z tego powodu, dla każdej wymienionej klasy modeli wyznaczone zostały prognozy dla kilku wersji modelu. Literatura zawiera opracowania przedstawiające różne podejścia do dni specjalnych w mo- delowaniu popytu na energię elektryczną. Weron8 omówił sposoby korygowania outlierów oraz uzupełniania braków danych. Harvey i Koopman w aplikacji swojej metody posłu- żyli się danymi ze skorygowanymi obserwacjami dla dni specjalnych9. Dni świąteczne i święta państwowe zostały przez nich zastąpione średnią z odpowiednich dni tygodnia po- przedniego i następnego. Ramanathan i inni10 nie opisali żadnych przekształceń swoich danych. Zdecydowałem oszacować każdy z modeli w eksperymencie był szacowany dla trzech wa- riantów bazy danych. Pierwsza z nich zawierała oryginalne i w żaden sposób nieskorygo- wane dane. W drugim zbiorze danych obserwacje odstające zostały zastąpione uzyskanymi dla nich wartościami dopasowanymi. Ze względu na to, że analizowane dane są bardzo mocno zagregowane w przekroju terytorialnym - łączne zapotrzebowanie dla całego kraju - analiza graficzna nie wskazała potencjalnie odstających obserwacji. W przypadku Polski wyłączenia większych obszarów i pozbawienie większej liczby odbiorców dostępu do elek- tryczności związane jest zwykle z silnymi wiatrami. Za potencjalnie mające największy wpływ i widoczne być może widoczne w skali makro są okresy kilku dni następujących po przejściu orkanu. W okresie 2002-2014 w Polsce odnotowano 11 przypadków tych silnych wiatrów, a informacje o terminie uderzenia orkanów przedstawione są w tablicy 5.1. Za okres potencjalnie widocznego wpływu orkanu na popyt i podaż energii uznane zostały ar- bitralnie 3 pełne dni po zjawisku. W trzecim typie zbioru danych dni specjalne, związane ze świętami kościelnymi i państwowymi, zostały zastąpione wartościami dopasowanymi z 8Weron R., 2006, Modeling and Forecasting..., str. 69. 9Harvey A.C., Koopman S.J., 1993, Forecasting Hourly Electricity..., str. 1232. 10Ramanathan R., Engle R., Granger C.J.V., Vahid-Araghi F., Brace C., 1997, Short-run Forecast... 5.2. METODA DNIA PODOBNEGO 95 na 1 okres, w tym przypadku tydzień, metoda naiwna pozwala wyznaczyć wartość prze- widywaną, która w stosunku do wartości rzeczywistej różni się o 0, 456%. Spośród trzech analizowanych metod najlepsze wyniki przy sporządzaniu prognoz na tydzień w przód uzy- skać można z metody naiwnej. Odmienne wnioski zostały otrzymane dla przewidywania wartości zapotrzebowania na energię na 2 i 3 tygodnie naprzód, dla których najlepszym sposobem okazał się wybór najbardziej podobnego dnia w 6-miesięcznej historii obserwa- cji. Jakością dopasowania w okresie czasu letniego i zimowego jest zgodna z wynikami uzyskanymi dla całej próby. Zaznaczyć jednak należy, że metoda dnia podobnego pozwala dokładniej szacować oczekiwaną wielkość zjawiskaw okresie obowiązywania czasu letniego. Uwzględnienie informacji o orkanach w Polsce nie pozwala na istotne poprawienie jakości dopasowania, gdyż wartości MAPE są na niemalże identycznym poziomie. Wydaje się, że najrozsądniejszym sposobem pominięcia wpływu dni specjalnych na pro- gnozy dla metody dnia podobnego jest usunięcie obserwacji za te dni z bazy danych. Dla najprostszej wersji procedury błąd MAPE spada o 0, 1 punktu procentowego dla pro- gnoz jednookresowych i praktycznie nie zmienia się dla pozostałych horyzontów. Większe spadki obserwować można dla dwóch pozostałych wariantów, ponieważ były one równe 1 3 i prawie 1 2 wartości błędów dla oryginalnych danych. Na podstawie uzyskanych wyników należałoby wyciągnąć wniosek, że w przypadku niedysponowania modelem ekonometrycz- nym najlepszą procedurą sporządzania prognoz byłoby wykorzystanie metody naiwnej do wyznaczania prognoz na 1 okres i metod dnia podobnego z przeszukiwaniem 6 poprze- dzających miesięcy przy prognozach dwu- i trzyokresowych. Pozornie zaskakujący wniosek o wyższości wyboru dnia najbardziej podobnego do progno- zowanego spośród 6 poprzednich miesięcy może być wytłumaczony konkluzjami z rozdzia- łów 2 i 3. Obserwowana tendencja wzrostowa w zapotrzebowaniu na energię elektryczną i widoczny wpływ stanu gospodarki tłumaczą powyższy rezultat. Selekcja historycznej obserwacji o zbliżonych charakterystykach na podstawie całej próby może teoretycznie do- prowadzić do wybrania obserwacji sprzed kilkunastu lat i z okresu o odwrotnym stanem koniunktury. Oznacza to, że jednym z kryteriów doboru najbardziej podobnej obserwacji powinno być podobieństwo fazy cyklu koniunkturalnego. 96 ROZDZIAŁ 5. PROGNOZOWANIE POPYTU NA ENERGIĘ ELEKTRYCZNĄ Ta bl ic a 5. 2: B łę dy M A P E dl a ró żn yc h w er sj im et od y dn ia po do bn eg o w % O gó łe m C za s zi m ow y C za s le tn i (1 dz ie ń) (2 dz ie ń) (3 dz ie ń) (1 dz ie ń) (2 dz ie ń) (3 dz ie ń) (1 dz ie ń) (2 dz ie ń) (3 dz ie ń) D an e or yg in al ne M et od a na iw na 0, 45 6 0, 87 5 1, 16 7 0, 53 2 0, 88 6 1, 15 6 0, 38 3 0, 86 4 1, 17 9 D zi eń po do bn y 6M 0, 61 0 0, 67 9 0, 61 2 0, 75 1 0, 76 0 0, 75 0 0, 47 3 0, 46 9 0, 47 6 D zi eń po do bn y 1, 31 0 1, 30 9 1, 37 6 1, 45 0 1, 44 6 1, 50 8 1, 17 4 1, 17 4 1, 24 6 D an e w yg ła dz on e M et od a na iw na 0, 45 9 0, 87 5 1, 16 6 0, 53 9 0, 88 7 1, 15 2 0, 38 3 0, 86 4 1, 17 9 D zi eń po do bn y 6M 0, 61 0 0, 61 1 0, 61 1 0, 75 0 0, 75 7 0, 74 7 0, 47 3 0, 46 9 0, 47 6 D zi eń po do bn y 0, 10 5 0, 11 1 1, 12 1 1, 15 3 1, 20 9 1, 28 3 0, 94 2 1, 00 9 0, 96 3 D an e be z dn i sp ec ja ln yc h M et od a na iw na 0, 35 8 0, 81 5 1, 12 1 0, 36 9 0, 79 2 1, 07 8 0, 34 7 0, 83 5 1, 16 2 D zi eń po do bn y 6M 0, 43 4 0, 44 5 0, 45 3 0, 44 4 0, 46 6 0, 47 7 0, 42 5 0, 42 5 0, 42 8 D zi eń po do bn y 0, 72 0 0, 75 0 0, 74 7 0, 72 7 0, 74 6 0, 76 9 0, 71 4 0, 75 5 0, 72 5 H or yz on t pr og no zy na 1, 2, 3 na st ęp ny dz ie ń oz na cz on y je st ja ko (1 dz ie ń) ,( 2 dz ie ń) ,( 3 dz ie ń) . Źr ód ło : O pr ac ow an ie w ła sn e. 5.3. MODEL REGRESJI 97 5.3 Model regresji Pierwszym wyborem w modelowaniu, oczywiście, jeśli tylko może być zastosowany, jest standardowy model regresji. Ramanathan i inni15 stwierdzili, że naturalną jest chęć wy- estymowania jednego, skomplikowanego modelu dla całego szeregu czasowego. Model szacowany dla wszystkich obserwacji ma postać zbliżoną do modelu (4.13)-(4.14), z nie- znaczną korektą dla funkcji określającej nieliniowy wpływ temperatury na zapotrzebowa- nie na energię, i przyjmuje postać yt = β0 + β1t+ β2godzinat + β3dzient + β4miesiact + β5godzinat ∗ dzient + β6t ∗miesiact + β7lowcloctat + β8sunt + β9visibilityt + β10dni_specjalnet + β11temp_medt + β121(temp_medt>16) ∗ tempt + εt. (5.1) Powyższa forma funkcyjna odpowiada modelowaniu funkcją ciągłą i przedziałami liniową. Test Ramsey’a, nazywany też testem Reset, służy do badania poprawności formy funk- cyjnej w modelu regresji. W teście Ramsey’a przeprowadza się regresję reszt et na zbiorze potęg wartości dopasowanych z modelu postaci yt = xtβ + α1ŷ 2 t + . . .+ αpŷ p+1 t + ηt. (5.2) Weryfikuje się następnie hipotezę łączną H0 : α1 = . . . = αp = 0 za pomocą statystyki NR2 −→ χ2 p. Statystyka testowa uzyskana dla pełnego modelu oszacowanego na próbie od 1 stycznia 2002 r. do 10 października 2013 r. wyniosła 2608,75, przy 5% wartości krytycznej 2,605. Należy zatem odrzucić hipotezę zerową o poprawności formy funkcyjnej modelu. Jeden model dla wszystkich kwadransów, oprócz negatywnej weryfikacji testem Ramsey’a, nie spełnia wymogu poprawności formy funkcyjnej sprawdzanego także przy zastosowaniu testu Chowa. Test ten pozwala weryfikować hipotezę o stabilności parame- trów w podpróbkach. Konstrukcja testu pozwala analizować równość oszacowań w wielu podpróbkach. Test Chowa do badania stabilności parametrów w wielu podpróbach nie jest zaprogramowany w programie Stata. Odpowiednie wyliczenia zostały wykonane samo- dzielnie. Do przeprowadzenia weryfikacji hipotezy o równości oszacowań dla wszystkich kwadransów niezbędne jest otrzymanie sumy kwadratów reszt dla modelu pełnego (SR) oraz dla modeli oszacowanych dla poszczególnych podprób, których suma oznaczana jest jako S. Statystyka testowa ma, w tym przypadku, postać F = (SR − S)/[K(m− 1)] S/(N −mK) ∼ FK(m−1),N−mK , (5.3) gdzie K jest liczbą zmiennych w modelu (K = 130), N = 412845 - liczbą obserwacji, natomiast m - liczbą analizowanych podpróbek. Uzyskana wartość statystyki testowej F 15Ramanathan R., Engle R., Granger C.J.V., Vahid-Araghi F., Brace C., 1997, Short-run Forecast..., str. 165