Rozdział 6 Bayesowska analiza statystyczna w badaniach koniunktury gospodarczej

Bayesowska analiza statystyczna w badaniach koniunktury gospodarczej

W niniejszym rozdziale opiszemy wnioskowanie statystyczne oparte na zasadzie Bayesa. Praktyczne stosowanie metod bayesowskich nie jest niczym nowym, a jest powszechniejsze tym bardziej, im ściślejszy jest związek statystycznej teorii i praktyki.

Niewątpliwie natomiast mniej niż metody klasyczne nadają się one do samodzielnego i rutynowego stosowania przez niespecjalistów. Elementarne podstawy statystycznej analizy bayesowskiej w zastosowaniu do badań koniunktury zostały opisane w pracy Męczarskiego (1998), a ogólnie na przykład w podręczniku Jóźwiak i Podgórskiego (1999), zaś w bardziej zaawansowanej postaci w monografii DeGroota (1983) i licznych podręcznikach statystyki matematycznej (Bartoszewicz, 1989).

6.1. Podstawy

Dane statystyczne gromadzone w badaniach ankietowych rozpisywanych w celu analizy koniunktury gospodarczej są w bezpośredni i naturalny sposób modelowane za pomocą rozkładu wielomianowego o trzech klasach M (p x, p 2, P i). Jest to rozkład typu dyskretnego o funkcji prawdopodobieństwa

P (X j — x x, X 2 — x2, X 3 — x3) — n\

*l!*2!x3!P \ P l P l > ^(6.1) gdzie xx+ x 2 + x3 = n, p x+ p 2 + p 3 =\, p, >0, i = 1,2,3. Liczby P \,P2->Pj, t0 odpowiednio frakcje (interpretowane też jako prawdopodobieństwa) odpowiedzi

„wzrost”, „bez zmian” i „spadek” , a *,,jc2,jt3 to liczby uzyskanych w badaniach odpowiedzi na pytanie o dany aspekt działalności gospodarczej badanych jednostek.

Wnioskowanie statystyczne polega na szacowaniu P \ , P i , p 2 w drodze estymacji punktowej lub przedziałowej oraz na weryfikacji hipotez statystycznych dotyczących tych parametrów oraz pewnych funkcji tych parametrów, dla potrzeb opisu, analizy zróżnicowania w czasie i przestrzeni oraz dla potrzeb prognozowania.

Przypominamy podstawowe parametry rozkładu M ( p x, p 2, p 3):

EX, - npi,

Z ) 2 X ; = 7 7 /7 , ( 1 - / ? , ) ,

Cov(X i, X ; ) = -7?p,P j, z * y, z, 7 = 1,2,3.

(6.2)

Klasyczna estymacja parametru p t metodą największej wiarogodności (ENW) i metodą momentów daje estymator nieobciążony o minimalnej wariancji (ENMW), o naturalnej postaci

6.2. Estymacja

Rozkład a priori odzwierciedla wiedzę a priori o parametrach rozkładu obserwacji.

Może również być wyrazem niewiedzy a priori i nieraz bywa wybierany w sposób formalny lub automatyczny, podobnie jak funkcja straty w zagadnieniach decyzyjnych.

Typowy wybór rozkładu a priori to tzw. rozkład sprzężony do rozkładu obserwacji.

Jest to cecha matematyczna (zob. np. DeGroot, 1983), dzięki której rozkład a posteriori ma tę samą postać analityczną, co rozkład a priori, tylko ze zmodyfikowanymi parametrami. Dobór parametrów rozkładu a priori omówimy w osobnym podrozdziale.

Rozkładem a priori parametrów P ^ Pi^Pt, sprzężonym do rozkładu M (/?], /?2, /?3) jest rozkład Dirichleta D{a{,a 1,a-i ). Jest to rozkład typu ciągłego o funkcji gęstości

gdzie A = (# ,,a2, a3) jest wektorem parametrów o dodatnich współrzędnych, a T oznacza funkcję gamma Eulera.

P i (6.3)

n Wtedy

(6.4)

r(<3j +^2 + «o)

P\ p2z p3

< 3 . - 1 - 1-1 -1 dla P]+ P j + P3 = h

f A P V P^Pri) =

oroz p. > 0, (6.5)

0 w przeciwnym przypadku,

Zauważmy, że parametry rozkładu M ( p l, p 2, p 3) mają sens prawdopodobieństw.

Odpowiada temu rozkład Dirichleta, ponieważ jego gęstość jest różna od zera tylko wtedy, gdy

P\ +P2 + Pi = Pi 2 0 ,i = 1,2,3,

a zatem wtedy, gdy sens prawdopodobieństw jest zachowany. Takiej własności należy oczekiwać od każdego rozkładu a priori dla parametrów rozkładu M ( p x, p 2, p 3) .

Rozkład a posteriori jest także rozkładem Dirichleta, a mianowicie rozkładem D(a] + x l,a-) + x2,a 3 + x3). Przy typowych założeniach statystycznego zagadnienia decyzyjnego z kwadratową funkcją straty estymator optymalny (bayesowski) jest to wartość oczekiwana w rozkładzie a posteriori. Ogólnie dla rozkładu Dirichleta

D{cx, c2,c3) mamy

EYi = --- ---, / = 1,2,3, c, + c2 + C3

więc estymator bayesowski parametru /?, ma postać

« a, + x, ci: a X: n . .

p ” = —— - = —--- h —---, gdzie a = a{ + a2 + a2. (6.6) a + n a a + n n a + n

Widać, że jest to średnia ważona wartości oczekiwanej a priori i estymatora o największej wiarogodności.

W zależności od wyboru parametrów rozkładu a priori uzyskujemy estymatory o nieco różnych postaciach (zob. Męczarski, 1998). Rozkład Z)(l,l,l) jest rozkładem równomiernym (jednostajnym) na zbiorze

S={(P\ >P2i Pl )Z R 3 : Pl + Pl +P3 = h P \ ’P l ’P3 ^°}>

który jest trójkątem na płaszczyźnie o równaniu p\ + p 2 + p 3 = 1 w przestrzeni R ’ leżącym w I oktancie układu współrzędnych. Rozkład równomierny jest przyjętym, choć uproszczonym sposobem wyrażania braku informacji a priori w postaci przekonania o „równouprawnieniu” wszystkich wartości parametru ze zbioru S . Głębsze rozważania (zob. Box, Tiao, 1973) doprowadzają do rozkładu jako najwłaściwszego z pewnego punktu widzenia na odzwierciedlanie „ignorancji”

a priori. Jest to tak zwany rozkład wzorcowy (ang. reference prior) dla rozważanego problemu (pojęcie to jest definiowane z wykorzystaniem rozważań o zasobie informacji). Na koniec uzyskany na gruncie ogólnej teorii estymator minimaksowy (minimalizujący maksimum ryzyka), a więc liczący się z najmniej pomyślnym doborem rozkładu a priori, okazuje się estymatorem bayesowskim przy rozkładzie

a priori Jest to tak zwany najmniej korzystny rozkład a priori, czyli taki rozkład a priori, że odpowiadający mu estymator bayesowski ma największe ryzyko spośród wszystkich estymatorów bayesowskich. Zatem z przyczyn teoretycznych rozpatrujemy trzy ze sprzężonych rozkładów a priori spośród symetrycznych rozkładów Dirichleta D ( a ,a ,a ): jednostajny, wzorcowy i najmniej korzystny. Gdy badamy wpływ zmienności parametru a e ( { ’ 3 ) na wartość estymatora, to okazuje się, że nie jest on znaczny. Estymator jest dość stabilny w warunkach zmienności rozkładu a priori w zakresie symetrycznych rozkładów Dirichleta. Ilustruje to poniższy przykład.

Przykład 6.1. Na pytanie o wielkość zatrudnienia w budownictwie w przekroju wszystkich przedsiębiorstw w I kwartale 2001 r. z reprezentatywnej próby n = 479 przedsiębiorstw otrzymano następujące odsetki odpowiedzi: wzrost: 3,7%, brak zmian:

35,7%, spadek: 60,6%. Wartości te, zgodnie ze wzorem (6.3), dają nam wartości ENMW, EN W i estymatora metodą momentów odpowiednich frakcji w całej populacji. Estymatory bayesowskie przyjmują wartości zależnie od rozkładu a priori:

dla rozkładu a priori wzorcowego p n = (0,037; 0,359; 0,604), dla rozkładu a priori jednostajnego p J = (0,037; 0,361; 0,602), dla rozkładu a priori najmniej korzystnego, estymator minimaksowy

p M =(0,036; 0,385; 0,579)).

Dla rozkładu najmniej korzystnego różnica między EMNW a estymatorem bayesowskim jest największa, ponieważ waga wartości oczekiwanej a priori równej

l we wzorze (6.6) jest w tym przypadku największa:

a

a + n 0,044.

yfn + n

Na rysunku 6.1 przedstawiamy porównanie bayesowskich ocen frakcji odpowiedzi

„wzrost” przy trzech rozkładach a priori w ciągu całego okresu badania.

Widoczne jest, że wybór rozkładu a priori spośród trzech rozpatrywanych nie ma dużego wpływu na wyniki estymacji. Jedynie estymator bayesowski na podstawie najmniej korzystnego rozkładu a priori różni się nieco wyraźniej od pozostałych.

6.3. Empiryczna estymacja bayesowska

Metoda empiryczna w estymacji bayesowskiej jest jednym ze sposobów pokonywania trudności z wytypowaniem odpowiedniego rozkładu prawdopodobieństwa jako rozkładu a priori. Dysponujemy najpierw rozkładem obserwacji o gęstości (lub funkcji prawdopodobieństwa w przypadku rozkładu typu dyskretnego) oznaczonej przez /( x |i 9 ) , traktując go jako rozkład warunkowy obserwacji lub statystyki X pod warunkiem znajomości & , oraz rozkładem a priori 0 gęstości 7t(3). Zgodnie z definicją rozkładu warunkowego konstruujemy łączny rozkład wektora losowego złożonego z obserwacji i parametru, który - jak wiadomo - w paradygmacie bayesowskim jest losowy, o postaci

p(x,&) = f ( x \ & y n ( 9 ) . (6.7) Według podstawowych reguł rachunku prawdopodobieństwa otrzymujemy gęstość rozkładu brzegowego obserwacji X

mn (x) = ^ p (x ,S )d 9 = ^f{x\& )jr(3)d& . (6.8)

© ©

Rozkład ten zależy od rozkładu a priori n : zarówno od jego postaci funkcyjnej, jak 1 od jego parametrów. Oznacza to w szczególności, że parametry rozkładu a priori (tak zwane hiperparametry) można oszacować na podstawie próby losowej. Przyjmujemy wówczas, że obserwacje mają rozkład (6.8) i takie oszacowanie wykorzystuje się następnie we wzorach opisujących estymatory bayesowskie. Podana metoda jest jasna intuicyjnie, ale budzi zastrzeżenia, ponieważ tak otrzymany rozkład nie jest rozkładem a priori w prawidłowym sensie, gdyż zależy od próby i nie spełnia wzoru Bayesa.

Otrzymujemy więc pewien estymator „niby-bayesowski”, który jest tylko opartym na analogii oszacowaniem podobnym do bayesowskiego i nie jest pewne, jakie ma własności. Prawidłowa empiryczna estymacja bayesowska to w istocie rodzaj prognozy z elementami wnioskowania nieparametrycznego (Robbins, 1955). Jeżeli jednak estymacji hiperparametrów dokonuje się na podstawie innych danych, niż estymacji bayesowskiej parametru *9, to poprawność jest zachowana (nie jest to z kolei „empiryzm bayesowski” w sensie Robbinsa). Problemy te znalazły precyzyjne i jasne omówienie w monografii Roberta (2001).

Rozważmy empiryczną estymację bayesowską dla danych trójmianowych z rozkładem Dirichleta jako rozkładem a priori. Rozkładem brzegowym obserwacji jest tak zwany rozkład wielomianowy-Dirichleta MD(tf,,<3 9,<33) (zob. Bernardo i Smith, 1994). Jest to rozkład wielowymiarowy typu dyskretnego o następującej funkcji prawdopodobieństwa w przypadku wymiaru 3:

P (X ^ x j , *2 »^ ^ ^3) — n\

M t a ] t a ]

/ 1 n ^4- n

l

* , t a t a ! „[«] ,[«]

A M

3

(6.9)

gdzie a i >0, + x 2 + x 3 = n, x , = 0,1,2,..., i = 1,2,3, przy czym 6[v] = ^ + l)>...-(Z? + 5 - l ) = r-^ - ^ .

F(b)

Postać ta jest dość skomplikowana, niemniej dla naszych obliczeń wystarczą momenty tego rozkładu. Niech

a.

' r - t wtedy

EX l =np„ D 2X, = C o v ( X „ X / ) = — --- P)P

1 + a \ + a

n + a

Próbując wykorzystać we wspomniany heurystyczny sposób metodę momentów, na podstawie porównania wartości oczekiwanych otrzymujemy równania

a x

— = - - , / = 1,2,3, gdzie a - a] + a7 + a3,

a n

a stąd oszacowania w postaci

<*,- =fl3 — , / = 1,2,3-

Empiryczny estymator „niby-bayesowski” przybiera wówczas postać

^,“ ' = — - ■ = " . ' = ,.2,3. (6.10) d + n n

Zatem otrzymaliśmy estymator o znanej w teorii klasycznej postaci i własnościach, nie będący jednak estymatorem bayesowskim względem ustalonego rozkładu a priori.

Szacując prawidłowo rozkład brzegowy (6.8) - (6.9) na podstawie danych z przeszłości napotykamy trudności: potrzebujemy oszacowania wariancji na podstawie przynajmniej dwuelementowej próby z ustalonego rozkładu wielomianowego, a tym w omawianych badaniach nie dysponujemy, ponieważ każdy wynik opiera się na innej liczbie odpowiedzi n. Wykorzystamy zatem oszacowanie typu empirycznego bezpośrednio dla rozkładu a priori, w ślad za Bergerem i Chenem (1993). Chodzi o to, że wartości oczekiwane a priori i brzegowe są identyczne w sytuacji rozkładu wielomianowego obserwacji i sprzężonego rozkładu a priori.

Przyrównujemy mianowicie oszacowania />, , i = 1,2,3, do składowych wektora wartości oczekiwanej rozkładu a priori D(ax, a2, a 3), a oszacowania wariancji

P i , i = 1,2,3, do odpowiednich elementów macierzy kowariancji tego rozkładu.

Wtedy:

skąd otrzymujemy

a _ 2 a ,( a - a , ) Pi ’ Si 2 / n ’

a a (a + 1)

_ P\ _ Pi

a\ ~ a 3 . » a2 ~ a3 . » a3

Pl P

/ = 1,2,3,

= { n - \ ) p :„

(6.11)

czyli

at = (n - 1 ) Pj, i = 1,2,3, skąd a = n - 1,

przy czym «jest liczbą odpowiedzi z okresu, z którego pochodzą oszacowania.

Aby przedstawione rozumowanie i wyniki były poprawne, oszacowania p, służące do ustalenia wartości aj nie mogą być obliczone na podstawie tych danych, które wykorzystujemy do wnioskowania o parametrach p ;, mając już wartości hiperparametrów. Wydaje się, że oszacowania warto konstruować, korzystając z danych z przeszłości (sprzed roku, wobec sezonowości) lub z prognoz respondentów ankiety na rozpatrywany kwartał. Poniższy przykład ilustruje wielkość różnic między estymatorami opartymi na danych z przeszłości a opartymi na prognozie.

Przykład 6.2. Ponownie rozpatrujemy wielkość zatrudnienia w budownictwie w przekroju wszystkich przedsiębiorstw. W I kwartale 2000 roku uzyskano

następujące frakcje odpowiedzi „wzrost”, „bez zmian” i „spadek”: 4,8%, 40,8%

i 54,4% (odpowiednio). Analogiczne prognozy na I kwartał 2001 roku podawane w IV kwartale 2000 dały wyniki: 1,8%, 39,9% i 58,3%. Zatem estymator bayesowski stanu z I kwartału roku 2001 dla rozkładu a priori o parametrach obliczonych na podstawie stanu sprzed roku ma wartości

P przesz = (0,043; 0,384 ; 0,573),

a na podstawie prognozy z poprzedniego kwartału p prog = (0,028; 0,377 ;0,595).

Na rysunku 6.2 przedstawiamy wartości estymatorów frakcji odpowiedzi „wzrost”

dla okresu badania od 1 kwartału 1995 roku (rok po rozpoczęciu badań). Można zauważyć, że w przypadku estymatorów bayesowskich typu empirycznego różnica między wartością takiego estymatora a estymatora NW odzwierciedla różnicę między stanem obecnym a przeszłością lub prognozą.

Rys. 6.2. Estymatory empiryczne bayesowskie oparte na prognozie oraz na stanie sprzed roku

6.4. Rozkłady predykcyjne - predykcja bayesowska

Swoisty schemat wnioskowania bayesowskiego, jakim jest wnioskowanie predykcyjne, polega na potraktowaniu rozkładu a posteriori | *) jako rozkładu apriori. Wartości obserwacji X stają się dodatkowymi parametrami, natomiast wnioskowanie dotyczy rozkładu przyszłej obserwacji Z . Mając jej rozkład warunkowy / (z | $), budujemy dla niej rozkład brzegowy w następujący sposób:

7r(z I x) = f f ( z I S)k{3 I x)d&

,

©

Zatem rozkład predykcyjny (predyktywny) Z jest to rozkład warunkowy tej zmiennej losowej pod warunkiem, że zaobserwowano X. Trzeba wspomnieć, że w lite­

raturze przedmiotu określa się czasem jako rozkład predykcyjny (predictive distributioń) zwykły rozkład brzegowy obserwacji mK(x) określony wzorem (6.8).

Omawiany teraz rozkład nazywa się wtedy rozkładem predykcyjnym a posteriori (posterior predictive distributioń).

Jeżeli rozkład obserwacji X = ( X x, X 2, X 3) jest wielomianowy M n( p x, p 2, p 2), rozkład obserwacji prognozowanej Z = (Z ,,Z 2,Z 3) jest rozkładem wielomianowym M r(Pi,p->,p3)-> a rozkład a priori rozkładem Dirichleta D(ax, a2, a3), to rozkład predykcyjny obserwacji Z jest typu wielomianowy-Dirichleta

MD(ax + x x, a2 + x 2,a 3 + jc3), którego wartość oczekiwana daje prognozę

o znanej postaci, identycznej ze zwykłym estymatorem bayesowskim przy rozkładzie a priori Dirichleta. Zasadą prognozy tego typu jest przyjęcie rozkładu a posteriori obliczonego na podstawie przeszłych obserwacji jako rozkładu a priori dla parametru rozkładu wielkości prognozowanej, stąd postać prognozy. Zaskakująca prostota jest też oczywiście skutkiem braku założeń o zależności w czasie w modelu statystycznym.

6.5. Bayesowska weryfikacja hipotez

6.5.1. Uwagi ogólne

Konstrukcja testów bayesowskich różni się od znanej dla testów budowanych na podstawie lematu Neymana-Pearsona. Tak samo definiujemy hipotezy statystyczne:

natomiast istnieją różne sposoby dalszego postępowania. Może ono prowadzić do konstrukcji zbioru krytycznego podobnie jak dla testów klasycznych. Mianowicie, definiujemy funkcję straty L{9,d) związanej z decyzją d o akceptacji i odrzucaniu poszczególnych hipotez; na przykład dla decyzji d$ „przyjąć hipotezę zerową” oraz

d x „przyjąć hipotezę alternatywną”

(6.13) a + n

H 0 : 9 e ®0, H X : 9 e ® x,@0 u ® x = 0 , © on 0 , = 0 , (6.14)

, Z (iM ,) =

O dla # e 0 , c, dla d e 0 O

Test bayesowski, czyli optymalny, powinien minimalizować ryzyko a posteriori Rn ( d , X ) decyzji d . Odrzucenie hipotezy zerowej przez ten test następuje wtedy, gdy

gdzie 0 o |x) jest prawdopodobieństwem a posteriori prawdziwości hipotezy zerowej, a więc wielkością obliczaną na podstawie obserwacji.

Procedura testowania może jednak polegać na innego rodzaju analizie prawdopodobieństw a priori i a posteriori prawdziwości obu hipotez. Wersja najbardziej uproszczona polega na podjęciu decyzji na podstawie wartości prawdopodobieństwa a posteriori hipotezy H 0. Jest to duże uproszczenie, a odrzucenie hipotezy zerowej tylko dlatego, że jej prawdopodobieństwo jest mniejsze od i , może być pomyłką. W istocie znaczenie ma fakt, czy prawdopodobieństwo hipotezy zerowej a posteriori jest większe niż a priori. Oznacza to wtedy, że dane przemawiają na korzyść hipotezy. Pełniejszej analizy dokonuje się za pomocą tak zwanego ilorazu szans. Iloraz szans a posteriori (posterior odds ratio) jest to wielkość

©i

przy czym dla rozkładów typu dyskretnego całkę rozumie się jako odpowiednią sumę.

Iloraz szans a priori jest to stosunek

Rx (dl, X ) < R x {d0, X) . Oznacza to, że obszar krytyczny ma postać

f f ( x | & )n (d )d d

= ^ ( t f G 0 o \x) = i

posl /^(z?G 0 , | x) J / ( x | d )n {d )d d P„( &e Q0 \x) _ £ {

Pn (t?£ 0 j ) {&)d&

©i Iloraz

nazywamy czynnikiem bayesowskim (.Bayes factor) na rzecz hipotezy H 0. Jest to też iloraz szans a posteriori, jeżeli / = 1, czyli gdy hipotezy są a priori jednakowo prawdopodobne. Odgrywa zatem rolę mnożnika zmieniającego iloraz szans a priori na podstawie obserwacji.

6.5.2. Weryfikacja hipotez o istotności różnic

Analizę wykorzystującą czynnik bayesowski i jego części składowe zastosujemy do weryfikacji hipotezy prostej H 0 :fl = i90 przeciw alternatywnej hipotezie złożonej

: d ■£ $0. Wtedy rozkład a priori musi być zadany tak, aby Pn {6 = d q) = Pq > 0, przy czym na zbiorze {# : d * $0} zadany jest rozkład a priori A , tak że

co należy rozumieć jako mieszankę rozkładu jednopunktowego skupionego w $ 0 oraz pewnego rozkładu prawdopodobieństwa na zbiorze W rozkładzie a posteriori określonym zgodnie z definicjąjako

mianownik jest funkcją gęstości rozkładu brzegowego obserwacji i ma postać

(6.15)

o

mn (x) = p 0f ( x | Ą ) + 0 “ Po

)J

I ^,

0

skąd wynika postać rozkładu a posteriori

(6.16)

Zauważmy, że prawdopodobieństwa a posteriori hipotez mają postać

P ( » = » \ x ) = p P * 0 | x) = - ...0 ...- - = 1

wn (jc)

- = l - P y (6.17) Zatem czynnik bayesowski określony jest jako

P\

l ~ P\

P_0 ' - P o

\ f ( x \ 9 ] A ( 9 ) d 9 0

/ ( * I « V

(6.18)

Jest to iloraz szans a posteriori na rzecz H 0 dla p 0 = j . Jest to też - przypominamy -mnożnik, którego trzeba użyć, aby obliczyć iloraz szans a posteriori przy danym ilorazie szans a priori. Zauważmy, że możemy posłużyć się czynnikiem bayesowskim pomijając sprecyzowanie prawdopodobieństwa p 0.

Przykład 6.3. W pracy Męczarskiego (1998) weryfikowano hipotezę dla prawdopodobieństwa p odpowiedzi „wzrost” na wybrane pytanie:

przyjmując rozkład a priori jako Pn (p = p Q) = /r0 oraz rozkład X na zbiorze (0,/?0) u ( / ?0,l) jako Beta{aA,a~,). Otrzymano prawdopodobieństwo a posteriori hipotezy zerowej w postaci

gdzie P(r,5') oznacza funkcję beta, n - liczbę wszystkich odpowiedzi, x - liczbę odpowiedzi „wzrost”. Dla przekroju określonego jako gałąź EKD 45.4 (prace wykończeniowe), gdzie n = 19, x = 6, weryfikowano hipotezę o braku istotnej różnicy między stanem z II kwartału 1994 r., p 0 = 0,421, a prognozą na III kwartał,

p = 0,316 . Wyniki były następujące.

1. Dla 7r0 = 3-, dla wartości parametrów rozkładu a priori beta a ] = a 2 = 1 (rozkład równomierny na przedziale (0,1)), iloraz szans a priori wyniósł / = ' , zaś iloraz szans a posteriori I pos, = « 1,242, oraz BFH « 2,484. Zatem brak było silnych podstaw do odrzucenia hipotezy zerowej;

2. Dla 7r0 = 5 , dla rozkładu a priori równomiernego na przedziale (0,1), iloraz szans a priori wyniósł I pr =1, iloraz szans a posteriori / , = «2,484, oraz BFHf « 2,484. Zatem nie było podstaw do odrzucenia hipotezy zerowej. Oczywiście w obydwu przypadkach czynnik bayesowski miał tę samą wartość;

Hq'-P = Pq przeciw H ] : p ^ p Q,

tt0 J 3 ( a x, a 2) p x0 0 - p 0 ) n x j o f i { ax + x , a ^ ■»■><*■ 22 + n - x )

3. Ponownie dla 71q = \ , ale dla wartości parametrów rozkładu beta a, =<x- _ \_n 2 ~ 2

(najmniej korzystny rozkład a priori), iloraz szans a priori miał wartość 7^. = ł,, iloraz szans a posteriori I = =0,957 oraz BFH =1,914, a więc czynnik bayesowski był większy od 1. Zatem mamy podstawę do nieodrzucenia hipotezy zerowej;

4. Dla 7T0 = \ , dla najmniej korzystnego rozkładu a priori, iloraz szans a priori

wyżej - różnica pochodzi od błędów zaokrągleń). Zatem nie było wyraźnych podstaw do odrzucenia hipotezy zerowej.

Testy bayesowskie są przydatne w przypadku przekrojów o małej liczebności, gdy nie da się wykorzystywać asymptotycznych rozkładów normalnych. W podanym przykładzie test okazał się dość „tolerancyjny”, aczkolwiek w kategoriach czynnika bayesowskiego przewaga hipotezy o równości nie była przygniatająca - 2,5 lub 2 do 1.

Z kolei dla większych liczebności test taki okazuje się zbyt restrykcyjny (zob.

podrozdział 6.10.2).

6.6. Wnioskowanie przy niesprzężonych rozkładach a p r i o r i

Rozkład asymptotyczny estymatora (/>,,/>3) przy n — czyli rozkład przybliżony dla dużych n, jest dwuwymiarowym rozkładem normalnym

N2((/?i,/>2)>■£)> gdzie macierz kowariancji ma postać

normalnym znacznie ułatwia obliczenia związane z testami, a zwłaszcza ze zbiorami ufności. Jeżeli jednak wprowadzimy rozkład a priori, to rozkład sprzężony jest również rozkładem normalnym. Parametry mające sens prawdopodobieństw mogłyby zatem przyjmować wartości spoza przedziału (0,1) i nie dawać w sumie 1.

wyniósł I pr = \ , iloraz szans a posteriori I t - ^ 5 = 1,915,a BFHq = 1,915 (jak

P \ 0 - P \ ) _ P\Pi

(6.19)

77 77

Saldo

B — P1 P i,

dla którego

EB = p ] - p 3 oraz D 2B = ■£ [/?, + p 3 - ( p ^ P s ) 2],

ma z kolei asymptotyczny rozkład normalny N ( E B , D 2 B) . Wartości parametrów p x i p 3 ponownie sprawiają opisany wcześniej kłopot.

Niezbędne ograniczenia można nałożyć w wygodny i naturalny sposób za pomocą rozkładu a priori, różnego od rozkładu normalnego. Zatem zamiast sprzężonego rozkładu a priori zastosujemy trójwymiarowy (w istocie dwuwymiarowy) rozkład Dirichleta D( ax, a - a x - a 3, a3) w przypadku salda oraz wektora ( P\ , Pi ) , ewentualnie rozkład beta dla pojedynczego parametru. Mianowicie, do wnioskowania o frakcji odpowiedzi „wzrost” p ] wykorzystujemy statystykę p x = o przybliżonym rozkładzie normalnym N ( p x, p ] (1 - /? ,) ) . Oznaczmy jego gęstość przez / {y | p ) , dla uproszczenia zapisu pomijamy tu wskaźnik przy p . Parametr /?, ma rozkład a priori beta o gęstości

r(tfi + an ) a - l a n -1

= (] ~ p ) 2 dla P G (°d)>

gdzie A = (ax, a2), a x, a2 > 0 . Rozkład a posteriori ma gęstość postaci

n A ( p \ y ) =

f ( y \ p ) n A ( p )

exp n ( y - p) 2 p ( \ - p )

a, - 1

1 ^{f ( y} I ^{p ) x A} (^{p ) dp}

0 - p)^{°2 ~}

(0,1)( P )

Jexpj-

0 | 2 p ( \ - p )

. (6.20)

o - \ a

\ P 1 ( I - / ? ) 2 dp

Do wnioskowania o saldzie, czyli różnicy między frakcją „wzrost” a frakcją

„spadek” wykorzystujemy statystykę B = p x - p 3 o rozkładzie asymptotycznym normalnym

(zob. Podgórska, 1998). Oznaczmy jego gęstość przez g(b\ p x, p 3) . Rozkład a priori dla wektora ( p x, p 3) jest rozkładem Dirichleta zapisanym w postaci

, . T ( a ,+ a 9 + a ,) -i-i,, sU na(P\>Pi ) = r T ~ m ~ ' x r / \ Pl r ( a x) r ( a2) r ( a3) V ~ P i ~ P ^ Pi

(P \, p 3) ma postać

s(b | p , p )7t (p p ) V y 3'

* A P ,’l \ \b) = ~A u l , y 3

l g(b \ Pv P3 ) n A ( p v P3 ) dP\dp3_{v j} 3j A ^ y y 3 J

\ J exp 0 0

exp- „ - P 3))2

2 p\ +pi ; (p\

(h — ( n - n

» v n

a

0 ~ P \ - P3) 2 P3 dP{dP3

Chodzi o rozkład a posteriori dla funkcji parametrycznej s = p x - p 3. Niech t = p 3 . Wtedy

Z powodu skomplikowanych całek wyniki estymacji punktowej i przedziałowej (zob. następny podrozdział) obliczane będą za pomocą metod numerycznych, także symulacyjnych, a zwłaszcza typowych dla zadań obliczeniowych statystycznej analizy bayesowskiej markowowskich metod Monte Carlo (Markov chain Monte Carlo).

6.7. Bayesowskie przedziały ufności

Do analizy salda zastosujemy estymację przedziałową. W bayesowskiej analizie statystycznej rozpatruje się dwa typy estymacji za pomocą zbiorów ufności:

bayesowskie „zbiory wiarogodne” (Bayesian credible sets; polska nazwa nie jest powszechnie ustalona) oraz zbiory o największej gęstości a posteriori (highest posterior density sets - HPD sets).

Przedziały wiarogodne są to przedziały oparte na kwantylach rozkładu a posteriori w następujący sposób:

o

gdzie P (# < 0* | x) = —, P_(jd > 6 X„ \ jc) = 1---- ,(6.22)

T

a więc krańce przedziału są odpowiednimi kwantylami rozkładu a posteriori, zaś poziom ufności przedziału, czyli jego prawdopodobieństwo a posteriori jest równe

1 - a .

Zbiory o największej gęstości a posteriori konstruuje się następująco. Wybiera się największą stałą c taką, aby zbiór

Wtedy W jest zbiorem o najmniejszej mierze (Lebesgue’a) na prostej rzeczywistej takim, że jego prawdopodobieństwo a posteriori jest równe y. Jeżeli rozkład a posteriori nie jest jednomodalny, to taki zbiór ufności nie musi być przedziałem.

Natomiast jeżeli rozkład a posteriori jest jednomodalny i symetryczny, to U = W . Z powodu złożoności obliczeń, a przede wszystkim braku rozwiązań danych za pomocą bezpośrednio użytecznego wzoru, zastosujemy metody numeryczne typu Monte Carlo.

Przykład 6.4. Rozpatrzmy portfel zamówień w przemyśle w przekroju wszystkich przedsiębiorstw. Na rysunku 6.3 przedstawiamy 95%-owe przedziały (pasy) ufności (czyli y = 0,95, a = 0,05 ) dla salda odpowiedzi w okresie od marca 1997 do września 1998. Jak widać, przedziały „wiarogodne” i typu HPD różnią się nieznacznie, co powinno wynikać z tego, że rozkład a posteriori jest dość regularny, przy wszystkich rozkładach Dirichleta wykorzystanych jako rozkłady a priori. Natomiast najkrótsze przedziały ufności otrzymuje się, wybierając a priori rozkład najmniej korzystny — ma on mniejszą wariancję (dla pojedynczej składowej) od rozkładu jednostajnego i wzorcowego. Można to zauważyć, porównując rysunki 6.3 i 6.4 (mimo, że drugi z nich wykonano dla dłuższego szeregu czasowego).

Wy = { 9 \ x { 9 \ x ) > c y} (6.23)

spełniał równość

(6.24)

Rys. 6.3. Pas ufności dla salda zamówień ogółem w przemyśle, obliczony na podstawie przedziałów „wiarogodnych” i przedziałów typu HPD; rozkład a priori najmniej korzystny

Rys. 6.4. Pas ufności dla salda zamówień ogółem w przemyśle, obliczony na podstawie przedziałów „wiarogodnych” i przedziałów typu HPD; rozkład a priori wzorcowy

6.8. Metody obliczeniowe. Markowowskie metody Monte Carlo

W tym podrozdziale omawiamy techniczną realizację obliczeń, aby zachować niezbędną samowystarczalność czytelniczą niniejszego rozdziału. Narzędzia obliczeniowe są bowiem niestandardowe, a polskojęzyczne źródła poruszające problemy markowowskich metod Monte Carlo rozproszone i ubogie (ostatnio wydano monografię Osiewalskiego, 2001). Wykorzystane źródła to książki Wieczorkowskiego i Zielińskiego (1997), Roberta i Caselli (1999) oraz Chena i in. (2000). Pierwsza

dotyczy tylko generacji liczb losowych, trzecia metod Monte Carlo dla potrzeb obliczeń w analizie bayesowskiej, a druga łączy całość zagadnień metod Monte Carlo.

Przeprowadzone przez nas obliczenia dotyczyły:

1° estymacji punktowej frakcji odpowiedzi „wzrost” przy rozkładzie a priori zarówno sprzężonym (w tym przypadku wystarczającym narzędziem był pakiet EXCEL), jak i normalnym, oraz konstrukcji bayesowskich przedziałów ufności,

2° bayesowskiej estymacji przedziałowej salda przy rozkładzie a priori normalnym.

6.8.1,Obliczenia dla frakcji odpowiedzi „wzrost”

We wnioskowaniu statystycznym o pojedynczej frakcji posługujemy się rozkładem dwumianowym obserwacji oraz rozkładem a priori beta. Do obliczeń przyjęto trzy różne rozkłady a priori beta z parametrami:

i \ n 1

a] = a2 = 1, ax = a2 = --- oraz a, = a2 = - .

2 2

Wartość oczekiwana w rozkładzie a posteriori parametru p , oznaczającego frakcję odpowiedzi „wzrost” o postaci

E * a( p \ y ) =

została obliczona numerycznie. Znalezienie „wiarogodnych” zbiorów ufności (Bayesicm credible set) oraz zbiorów o największej gęstości a posteriori wymagało zastosowania metod symulacyjnych. Wygenerowana próba M-elementowa z rozkładu a posteriori została uporządkowana w sposób rosnący. Elementy -ty oraz

7 -V

1 -ty są odpowiednio lewym i prawym końcem przedziału „wiarogodnego”

([#] oznacza część całkowitą liczby a ). Algorytm szukania przedziału o największej gęstości a posteriori (przedziału HPD) polega na znalezieniu wszystkich przedziałów o końcach j i j + [(1 - a ) M ] , gdzie j = 1,2,..., M - [(1 - a ) M ] i wybraniu spośród

nich tego o najmniejszej mierze (Lebesgue’a) na prostej. Można pokazać, że przy pewnych założeniach (Chen i in., 2000) procedura ta daje przedział zbieżny do przedziału HPD, zdefiniowanego we wzorze (6.23). W obliczeniach przyjmujemy a - 0,05.

Osobnym zagadnieniem jest generowanie próby z rozkładu a posteriori.

W przypadku parametrów ax = a2 = 1 oraz ax = a 2 = ^ - generowanie zostało przeprowadzone za pomocą podstawowej wersji metody eliminacji. Mianowicie, w celu otrzymania liczb losowych z rozkładu o gęstości / n a przedziale (a,b), ograniczonej przez pewną stałą d > 0 należy wygenerować dwie niezależne zmienne losowe: U\ z rozkładu jednostajnego U(a,b) i U 2 z rozkładu U(0, d). Jeżeli spełniony jest warunek U2< f ( U 1), to przyjmuje się X = U 1, a w przeciwnym przypadku procedurę powtarza się od początku (por. Wieczorkowski i Zieliński, 1997).

Zgodnie z zasadami metody eliminacji, dla parametrów ax = a2 = 1 generujemy próbę z rozkładu o gęstości / proporcjonalnej do

exp " ( y - p ) 2 p ( \ - p )

2 \

a dla parametrów ax - a2 = — proporcjonalnej do\ n

exp

>2 A

W - p ) ) 2 ■ n { y - p Y

2p{\ - p)

Funkcje te są ograniczone przez d = 1, poza tym p e (0,1), zatem przyjmujemy a = b — d — \ .

Dla parametrów ax = a2 = ]- generowanie odbywa się z rozkładu o gęstości proporcjonalnej do

V 2/?0 - p ) )

Jest to funkcja nieograniczona i należy zastosować wariant ogólny metody eliminacji. Wybiera się funkcję g taką, że z rozkładu o tej gęstości umiemy otrzymywać liczby losowe oraz spełniony jest warunek f ( x ) < c g ( x ) . Algorytm polega na generowaniu zmiennej X z rozkładu o gęstości g i zmiennej U z rozkładu

jednostajnego na przedziale (0,1) tak długo, aż zostanie spełniony warunek

• cUg(x) < f ( x ) i na przyjęciu jako wyniku generowania wartości X . Oczywiście

expf n ( y - p ) 2 ^ 1 _1

(/>(!-/>)) 2 < ( p ( \ - p ) ) 2 v 2P d - p)

zatem c - 1, a g-(jc)jest funkcją gęstości rozkładu Beta\

wygenerowanej wartości X ma postać U < exp

( a

( n( y - p ) 2 ^

2 p ( l - p )

6.8.2. Obliczenia dla salda

Obliczenia wartości salda zostały przeprowadzone za pomocą markowowskich metod Monte Carlo. Zastosowano algorytm zwany próbnikiem Gibbsa, który przez konstrukcję odpowiedniego łańcucha Markowa ^p , gdzie ^p ' = ( p [ , p 3 ), pozwala otrzymywać losowy wektor z rozkładu a p o s t e r i o r i (6.21) przy założeniu znajomości rozkładów brzegowych a p o s t e r i o r i parametrów /?, i ^{p 3}, to znaczy umiejętności generowania zmiennych losowych o takich rozkładach.

Algorytm próbnika Gibbsa, który podajemy tu w wersji zastosowanej przez nas, jest następujący (por. Robert, Casella, 1999). Zakładamy, że znane są (umiemy generować) rozkłady brzegowe n ( p x \ b , p ^ ) i 7 t { p ^ \ b , p x ) . Niech ( p Q\ , p l ) będzie zaobserwowaną wartością początkową. Generujemy ^{p x} z rozkładu o gęstości

7 t ( p x \ b i p l ) , uzyskując w ten sposób wartość p \; p 3 generujemy z rozkładu

o gęstości /r(/?3 | b , p \ ) , otrzymując wartość ^{p \ .} Kontynuując, ^{p \} otrzymujemy jako wynik generacji z rozkładu ^{x ( p x} \ ^{b , p \ )} itd. Po wykonaniu T iteracji mamy wektor

p l = ( p l , p l), będący realizacją łańcucha Markowa z prawdopodobieństwami

przejścia

P{p', P ' +X) = n (p y \ b , p i ' ) t r ( p 3

1

Wykazuje się, że przy 1 dążącym do nieskończoności ^{p '} dąży według rozkładu do losowego wektora { p x , p 3 ) o rozkładzie takim, jak rozkład a p o s t e r i o r i o gęstości

n ( p \ , p 3 \b). Do obliczeń przyjęto T = 1000. Należy zauważyć, że kolejne

losowej, procedura została powtórzona M = 1000 razy i na jej podstawie obliczono wartości oczekiwane parametrów p ] i p z oraz wartość salda jako ich różnicę.

Korzystając z próby otrzymanej dzięki zastosowaniu próbnika Gibbsa, wartości oczekiwane powyższych parametrów wyznaczono inną metodą, na podstawie rozkładów brzegowych. Przedział (0,1) został podzielony na dziesięć części, a następnie zliczano, ile elementów próby znalazło się w poszczególnych podzbiorach podziału. Na podstawie tak otrzymanego histogramu obliczano wartości oczekiwane parametrów p x i p 3 oraz salda, jako ich różnicę.

Podamy jeszcze szczegóły generacji przy różnych rozkładach a priori wraz ze wzorami.

Generacja z rozkładów brzegowych o gęstościach

była dokonywana dla rozkładu a priori Dirichleta o parametrach a, = a 2 =\ oraz otrzymywane p ' nie są niezależne. Ponieważ naszym celem jest uzyskanie próby

n A(P\ Ib , p 3) = n ( f t - ( P i - p 3))2 P\

2 P\ + P i ~ ( P \ ~ P:J2 \{ + P i ~ ( P \ ~ P i f

J(

1

\Uj-l n_ (fr~(/?i ~ P 3))2 ] ( 1

2 P\ + P3 ~(P\ ~ P 3 ) 2 jl

1 ( ^ " ( P l - P s ) ) 2 \ ( { _ P\ T

oraz

n A( P3\ b>P\) =

a] = a 2 za pomocą podstawowej wersji metody eliminacji, analogicznie jak dla obliczeń frakcji odpowiedzi „wzrost”.

Dla parametrów z/, = a 2 = - zastosowano, podobnie jak dla frakcji „wzrost”, wariant ogólny metody eliminacji. Funkcja gęstości rozkładu, z którego chcemy generować próbę, jest nieograniczona i proporcjonalna do

/ ( P i ) = exp - ⁿ

( 6 - ( p , - p 3 ) ) 2 1 Ti -

* j ^{\ ~ 2}

l

J

1_ftJ

gdzie f i ( \ , ' ) oznacza funkcję beta. Za funkcję g (/?,) przyjęto f f

g(P\) =

V l -•

1 -

vv p 3

P\ W

A P3 )) 1 -

P3 Ą \ \

dla p { e (0,1 - )

i jest to funkcja gęstości dla rozkładu j3( \ , ' ) na przedziale (0,1 - /?,).

Warunek akceptacji wygenerowanej wartości ma wówczas postać I / > e x p ( - " Jb-M-PA I

i 2 P \ + P l - ( P \ ~ P l ) )

W obliczeniach wykorzystano własne implementacje metod Monte Carlo, napisane w języku C oraz generator liczb losowych z rozkładu jednostajnego na przedziale (0,1) ULTRA Marsaglii i Zamana (1992) o okresie powyżej 10 ?ó (public domain dla zastosowań niekomercyjnych). Generator ten został opisany przez Wieczorkowskiego i Zielińskiego (1997).

6.9. Wyniki

Przedstawione w poprzednich częściach rozdziału metody zastosowaliśmy do skonstruowania estymatorów bayesowskich punktowych i przedziałowych. Badanie dotyczyło zmian wielkości zatrudnienia i portfela zamówień krajowych w budownic­

twie oraz zmian wielkości zatrudnienia i portfela zamówień ogółem w przemyśle.

Wielkość zamówień krajowych należy do najważniejszych bodźców rozwoju polskiego budownictwa, a wielkość zatrudnienia - do najważniejszych przejawów ożywienia gospodarczego w tej dziedzinie. Podobnie jest z wielkością zamówień ogółem w przemyśle, a i wielkość zatrudnienia odbija w znacznym stopniu jego tendencje rozwojowe. Podajemy wyniki w postaci wykresów z krzywą wartości

estymatorów bayesowskich i pasami ufności 95% obu typów. Opierają się one na asymptotycznych normalnych rozkładach obserwacji. Rezultaty obliczeń na podstawie dokładnego rozkładu wielomianowego odgrywają rolę uzupełniającą (brak tu zbiorów ufności). Metodami testów bayesowskich zbadano również trafność prognoz portfela zamówień krajowych w budownictwie.

Wszystkie wyniki zostały otrzymane na podstawie danych Instytutu Rozwoju Gospodarczego SGH.

6.9.1. Budownictwo

Na poniższych rysunkach przedstawiamy wyniki estymacji wybranych wielkości opisujących koniunkturę w budownictwie, występujących W legendzie wykresu używamy następujących skrótów: HPDI oznacza przedział o największej gęstości a posteriori, a CI - przedział „wiarogodny”, czyli oparty na kwantylach rozkładu a posteriori.

Rys. 6.5a. Zbiór ufności na poziomie 95% dla salda wielkości zatrudnienia w budownictwie; wszystkie przedsiębiorstwa; roczna agregacja danych

zatrudnienia w budownictwie; wszystkie przedsiębiorstwa; roczna agregacja danych

Rys. 6.6. Zbiór ufności na poziomie 95% dla salda wielkości zatrudnienia w budownictwie; wszystkie przedsiębiorstwa; dane nieagregowane

Rys. 6.7. Zbiór ufności na poziomie 95% dla salda wielkości zatrudnienia w budownictwie; sektor publiczny; dane nieagregowane

Rys. 6.8. Zbiór ufności na poziomie 95% dla salda wielkości zatrudnienia w budownictwie; sektor prywatny; dane nieagregowane

Kwartał: rok. miesiąc

199801 199804 199807 199810 199901 199904 199907 199910 200001 200004 200007 200010 200101 200104

— — — Dolny kraniec HPDI — — — Górny kraniec HPDI - - - Dolny kr. Cl - - - Górny kr. Cl Ocena sald

Rys. 6.9a. Zbiór ufności na poziomie 95% dla salda liczby zamówień krajowych w budownictwie; wszystkie przedsiębiorstwa; roczna agregacja danych

Kwartał: rok, m iesiąc

199801 199804 199807 199810 199901 199904 199907 199910 200001 200004 200007 200010 200101 200104

Ocena frakcji — — — Dolny kraniec HPDI — — — Górny kraniec HPDI - - - Dolny kraniec Cl - ■ • • • -Górny kraniec Cl

Rys. 6.9b. Zbiór ufności na poziomie 95% dla frakcji odpowiedzi „wzrost” o liczbie zamówień krajowych w budownictwie; wszystkie przedsiębiorstwa;

roczna agregacja danych

Kwartał: rok, miesiąc

199704 199707 199710 199801 199804 199807 199810 199901 199904 199907 199910 200001 200004 200007 200010 200101 200104

--- Dolny kraniec H P D I--- Górny kr. HPDI - — - - Dolny kr. Cl - - * - - Górny kr Cl Ocena salda

Rys. 6.10. Zbiór ufności na poziomie 95% dla salda liczby zamówień krajowych w budownictwie; wszystkie przedsiębiorstwa; dane nieagregowane

Kwartał: rok. miesiąc

199704 199707 199710 199801 1998C4 199807 199810 199901 199904 199907 199910 200001 200004 200007 200010 200101 200104

Rys. 6.11. Zbiór ufności na poziomie 95% dla salda liczby zamówień krajowych w budownictwie; sektor publiczny; dane nieagregowane

Kwartał: rok. miesiąc

199704 199707 199710 199801 199804 199807 199810 199901 199904 199907 199910 200001 200004 200007 200010 200101 200104

Rys. 6.12. Zbiór ufności na poziomie 95% dla salda liczby zamówień krajowych w budownictwie; sektor prywatny; dane nieagregowane

Rocznej agregacji danych dokonywano w następujący sposób. Niech Y/4 oznacza wielkość z okresu i po agregacji, a Yt - wielkość nieprzekształconą. Wtedy przyjmujemy

i-4+k

, k = 1,2,3,4, Yi = a]Yi_3 + + tf3yM + aĄYt , gdzieak = ---

3 + 2 + */-i +

przy czym oznacza liczbę odpowiedzi w okresie r (por. Podgórska, 1998).

Agregacja skraca zatem szereg czasowy o 3 kwartały. Dane o zatrudnieniu są dostępne od I kwartału 1994 roku, a o zamówieniach krajowych - od I kwartału 1997 (wcześniej były to zamówienia ogółem). Sytuacja pod względem zatrudnienia była najlepsza w 1997 roku. Na przełomie 1998 i 1999 r. w obu kategoriach nastąpił duży sezonowy spadek, który okazał się jak dotąd nie do odrobienia. Tendencja spadkowa dla zatrudnienia utrzymuje się, a w zasadzie utrzymuje się dla zamówień krajowych. Jest ona słabsza, jeżeli mierzyć j ą frakcją odpowiedzi „wzrost”. Oznaczałoby to, że szybciej pogarsza się sytuacja przedsiębiorstw średnio prosperujących (wcześniej udzielających odpowiedzi „brak zmian”). Przedsiębiorstwa publiczne szybciej redukowały zatrudnienie, a ich portfel zamówień krajowych przedstawia się nieco korzystniej niż dla sektora prywatnego.

6.9.2. Trafność prognoz w budownictwie

Zastosujemy zasady bayesowskiej weryfikacji hipotez opisane w podrozdziale 6.6, a w szczególności w przykładzie 6.3, do testowania zgodności prognozy z realizacją dla portfela zamówień krajowych w budownictwie. Dla każdego kwartału definiujemy frakcję prognoz „wzrost” jako p 0. Weryfikujemy hipotezę, że stan p z tego kwartału (frakcja odpowiedzi „wzrost”) jest równy p 0. Przeprowadzamy test dla kolejnych kwartałów, obliczając ilorazy szans i czynniki bayesowskie, a wyniki przedstawiamy na rysunku 6.13a.

199707 199710 199801 199804 199807 199810 199901 199904 199907 199910 200001 200004 200007 200010 200101 200104 Kwartał: rok, miesiąc

" ^ “ ■Czynnik bayesowski

Rys. 6.13a. Badanie trafności prognoz za pomocą testu bayesowskiego: wartości czynnika bayesowskiego na rzecz hipotezy o całkowitej dokładności prognoz

Okazuje się, że prognozy nie są trafne. Żadna wartość czynnika bayesowskiego na rzecz równości prognozy i realizacji nawet nie zbliżyła się do 1. Jednak w rzeczywisto­

ści dla dużych liczebności przekroju test tego typu jest zbyt precyzyjny (przykład 6.3 dotyczył obliczeń dla niezbyt licznych przekrojów). Porównajmy bowiem wartości prognoz z realizacjami, a następnie z wartościami estymatora bayesowskiego opartego na prognozie (Tablica 6.1 i rys. 6.13b). Prognozy rzeczywiście bywają niezbyt trafne, ale brak trafności nie osiąga wielkich rozmiarów, tym bardziej, że niektóre sezonowe zmiany są w zasadzie pewne. Test bayesowski w tym prostym przypadku jest więc poniekąd zbyt dobry, aby był przydatny. Można sądzić, że spełniłby lepiej swoje zadanie w subtelniejszych zadaniach weryfikacji modelu (por. 6.10.4).

Kwartał: rok. miesiąc

•Estymator bayesowski ■ ENMW/ENW A Prognoza

Rys. 6.13b. Badanie trafności prognoz: estymator klasyczny i bayesowski oraz wartości prognoz

Tablica 6.1. Wartości prognoz i estymatorów służące do konstrukcji rys. 6.13a-b Kwartał BF - stan

a progn. Prognoza Stan Est. b.

199707 0,00004 66,7 56,7 61,8962

199710 3E-07 42,7 56 49,29949

199801 0,002 25,9 33,7 29,67453

199804 0,355 15,6 16,9 16,2939

199807 0,0002 58,9 48,9 54,22637

199810 0,043 41,1 46 43,60729

199901 0,002 27,1 32,4 30,55024

199904 0,22 15,7 13,9 14,907

199907 1E-08 59,8 47,5 54,0374

199910 4E-07 40,9 52,5 46,8142

200001 0,001 20,4 26,5 . 23,5665

200004 0,013 16,1 12 14,0095

200007 2E-07 61,4 49,2 55,78142

200010 0,000007 39,6 49,9 44,92946

200101 0,004 19,8 25,4 22,53791

200104 0,015 11,5 7,6 9,516595

Źródło: obliczenia własne na podstawie danych IRG SGH.

6.9.3. Przemysł

Na poniższych rysunkach przedstawiamy wyniki bayesowskiej estymacji przedziałowej wybranych wielkości opisujących koniunkturę w przemyśle.

W celu zestawienia metod rysunki 6.14b i 6.17b ilustrują estymację punktową

na podstawie sprzężonych rozkładów

Rys. 6.14a. Zbiór ufności na poziomie 95% dla salda liczby zamówień w przemyśle;

wszystkie przedsiębiorstwa

Rys. 6.14b. Estymatory punktowe (NMW/NW i bayesowskie) dla salda liczby zamówień w przemyśle; wszystkie przedsiębiorstwa

Rys. 6.14c. Zbiór ufności na poziomie 95% dla frakcji odpowiedzi „wzrost” o liczbie zamówień w przemyśle; wszystkie przedsiębiorstwa

sektor publiczny

sektor prywatny

Rys. 6.17a. Zbiór ufności na poziomie 95% dla salda wielkości zatrudnienia w przemyśle; wszystkie przedsiębiorstwa

Rok, m iesiąc

-5 -10

£ -15 S -20 2 -25 W -30 -35 -40

.* ^ “ 1 ~I" 1 " _____ ł ~ * ł j r > v --- _ ---

: : : : : : :

---Dane — — —— Est. minimaksowy - - - Est. b. wg prognozy

12 okr. śr. ruch. (Dane) 12 okr. śr. ruch. (Est. m in im a k s o w y ) ---12 okr. śr. ruch. (Est. b. wg prognozy)

Rys. 6.17b. Estymatory punktowe (NMW/NW i bayesowskie) dla salda wielkości zatrudnienia w przemyśle; wszystkie przedsiębiorstwa

o wielkości zatrudnienia w przemyśle; wszystkie przedsiębiorstwa

Rok, miesiąc

oB' oS? ,SN (*N # oB' o'?’ (f? £t '- ^ ^ ^

£ N# & & & & ^ ^ r#

Rys. 6.18. Zbiór ufności na poziomie 95% dla salda wielkości zatrudnienia w przemyśle; sektor publiczny

i

i t

l {l

i

» i t t t

i

Rys. 6.19. Zbiór ufności na poziomie 95% dla salda wielkości zatrudnienia w przemyśle; sektor prywatny

Sytuacja przemysłu również odznacza się sezonowością, ale nie aż tak silną, jak w budownictwie. Nie dokonywaliśmy więc agregacji danych, poprzestając na wykresie trendu (ruchoma średnia) sporządzanym w pakiecie EXCEL. Dane gromadzone są od marca 1997.

Najlepsza sytuacja pod względem portfela zamówień panowała podczas jesiennego ożywienia w 1997 r. (tylko wtedy saldo było dodatnie). Zima i wczesną wiosną 1998/1999 (po kryzysie w Rosji) nastąpił bardzo głęboki spadek, który już się nie powtórzył, ale i nie został zrekompensowany. Sytuacja w pierwszej połowie roku 2001 jest gorsza niż w roku 2000. W terminach frakcji odpowiedzi „wzrost” otrzymujemy podobne wyniki. Ożywienie 1997 roku było większe w sektorze prywatnym, ostatnio zaś była lepsza sytuacja sektora publicznego.

Zatrudnienie jest stabilniejsze niż portfel zamówień, zwłaszcza w sektorze publicznym. Wyraźne są zimowe spadki na początku lat 1998 i 2001. W dłuższym okresie sytuacja wydaje się polepszać, choć w terminach frakcji odpowiedzi „wzrost”

obserwujemy tendencję spadkową (w sposób oscylacyjny). Zatem coraz więcej przedsiębiorstw podaje odpowiedź „bez zmian”, zwłaszcza spośród sygnalizujących spadek. W ostatnim okresie nieco lepiej przedstawia się zatrudnienie w sektorze prywatnym. Byłby to pomyślny objaw, ponieważ sektor prywatny jest mniej skłonny do zatrudniania zbędnych pracowników; można też sądzić, że nastąpiła już daleko posunięta racjonalizacja zatrudnienia, a dalsza oznaczałaby zwijanie działalności.

Zwróćmy jeszcze uwagę na rysunki 6.14b i 6.17b, z wykresami punktowych estymatorów bayesowskich skonstruowanych metodą empiryczną na podstawie prognoz. Porównanie ich z danymi, a ściśle mówiąc z ENW, daje wyobrażenie o różnicach między prognozą a realizacją. Okazuje się, że wzrost wielkości portfela

zamówień prognozowało zawsze wyraźnie więcej przedsiębiorstw, niż go później odnotowało - realizacja była zawsze gorsza od prognozy (nawet o kilkanaście punktów procentowych). Natomiast zatrudnienie zwiększało się częściej, niż było prognozowane. Jest to więc pewna możliwość oceny trafności prognoz, choć otrzymana jako efekt uboczny i nie tak sformalizowana, jak w podrozdziale 6.10.2.

6.9.4. Uwagi końcowe

W niniejszym rozdziale przedstawiliśmy zarówno metody wraz z przykładami, koncentrując się na własnościach statystyczno-matematycznych, jak i wyniki analiz ekonometrycznych dla danych z ankietowych badań koniunktury gospodarczej. Istnieją możliwości udoskonalenia opisanych technik przez wprowadzenie jeszcze nowych elementów (rozbudowa testów), a także poprawienia efektywności programów symulacyjnych przez polepszenie metod generowania liczb losowych (niektóre przebiegi symulacji Monte Carlo trwały dość długo). Zasadniczym krokiem naprzód, do wykonania w przyszłych badaniach, byłoby jednak zastosowanie metod bayesowskiej analizy procesów stochastycznych i ekonometrii do specyficznych szeregów czasowych pojawiających się w ankietowych badaniach koniunktury.

Wzorem mogłyby być procedury estymacji, testowania i predykcji rozwijane dla danych obserwowanych na rynkach finansowych, również i w Polsce, choćby w duchu nowej monografii Osiewalskiego (2001). Można przypuszczać, że bayesowska analiza regresji przydałaby się również do badania powiązań między odpowiedziami na poszczególne pytania ankietowe.

Literatura

[ 1 ] J. Bartoszewicz (1989), Wykłady ze statystyki matematycznej, PWN Warszawa.

[2] J. O. Berger, M.-H. Chen (1993), Predicting retirement pattern: prediction for a multinomial distributioń with constrained parameter space, The Statistician 42, 427-443.

[3] J. Bernardo, A. F. M. Smith (1994), Bayesian Theory, J. Wiley, New York.

[4] G. E. P. Box, G. C. Tiao (1973), Bayesian Inference in Statistical Analysis, Addison-Wesley, Reading, Massachusetts.

[5] M. H. Chen, Q. M. Shao, J. G. Ibrahim (2000), Monte Carlo Methods in Bayesian Computations, Springer-Verlag, New York.

[6] M. DeGroot (1981), Optymalne decyzje statystyczne, PWN Warszawa.

[7] J. Jóźwiak, J. Podgórski (1997), Statystyka od podstaw’ (wyd. IV), PWE Warszawa.

[8] G. Marsaglia, A. Zaman (1992), ULTRA 1.01, Florida State University, Tallahassee.

[9] M. Męczarski (1998), Zróżnicowanie koniunktury. Metody wnioskowania statystycznego, w: Statystyczne i ekonometryczne metody badania

ła-ótkookresowych zmian stanu gospodarki, IRG SGH, Warszawa; str. 31-63.

[10] J. Osiewalski (2001), Ekonometria bayesowska w zastosowaniach, Wyd. AE Kraków.

[11] M. Podgórska (1998), Statystyczna analiza wyników testu koniunktury, w:

Statystyczne i ekonometryczne metody badania krótkookresowych zmian stanu gospodarki, IRG SGH, Warszawa; str. 65-88.

[13] C. P. Robert (2001), Bayesian Choice (2nd ed.), Springer-Verlag, New York.

[14] C. P. Robert, J. Casella (1999), Statistical Monte Carlo Methods, Springer-Verlag, New York.

[15] H. Robbins (1955), An empirical Bayes approach to statistics, w: ThirdBerkeley Symposium Math. Statist. Prób., Vol. 1, 157-163.

[16] R. Wieczorkowski, R. Zieliński (1997), Komputerowe generatory liczb losowych, WN-T Warszawa.

Bayesian statistical analysis in investigations of business activity for survey data

A bstract

In Chapter 6 applications of Bayesian statistics to investigations of business situation are treated. One presents foundations of Bayesian statistical inference for the multinomial model, in particular of interval estimation and statistical tests, computational probleins of applications including Markov Chain Monte Carlo methods and sonie results for business indicators for construction and industry. The data eonie from the database of IRG SGH (Institute of Economic Development, Warsaw School of Economics).