Ekonometria repetytorium
1. Przedmiot i narzędzia ekonometrii
Ekonometria nauka o metodach badania ilościowych prawidłowości występujących w
zjawiskach ekonomicznych. Wykorzystuje do tego aparat rachunku prawdopodobieństwa i
statystyki matematycznej oraz algebrę liniową (rachunek macierzowy)
Nieznane czynniki losowe v
Badane
zjawiska
ekonomiczne
Zm.egzogeniczne X Zm.endogeniczne: Y=f(X)+e
Wpływ otoczenia na proces Wpływ procesu na otoczenie
opisują zmienne objaśniające opisują zmienne objaśniane
(egzogeniczne) X (endogeniczne) Y: Y=f(X)+e
e reprezentacja losowości opisu
Rys.1. Ekonometryczne ujęcie zjawisk ekonomicznych
Ekonometria zajmuje się poszukiwaniem zależności ekonometrycznych f(X) (tj.
deterministycznych powiązań ilościowych pomiędzy zmiennymi objaśniającymi i
objaśnianymi) oraz analizą probabilistyczną składowej losowej e zmiennych objaśnianych.
Literatura:
1. Henry Theil: Zasady ekonometrii, PWN, Warszawa, 1979
2. Zbigniew Pawłowski: Ekonometria, PWN, Warszawa 1969
3. Edward Nowak: Zarys metod ekonometrii zbiór zadań, PWN, Warszawa 1994
4. John Freund Podstawy nowoczesnej statystyki, PWE, Warszawa 1968
5. G.E.P.Box, G.M.Jenkins: Analiza szeregów czasowych, PWN, Warszawa, 1983
2. Podstawowe pojęcia rachunku prawdopodobieństwa
Literatura:
1. I.E. Brontsztejn, K.A.Siemeindiajew: Matematyka poradnik encyklopedyczny. Część
szósta Opracowanie danych doświadczalnych,. PWN, Warszawa 1986
2. Poradnik inżyniera Matematyka Rozdziały XXXII i XXXIII
3. J.Greń: Statystyka matematyczna modele i zadania, PWN, Warszawa 1982
Definicje intuicyjne: (Foralnie definicje można znależć np. w Poradniku 2)
Zdarzenie losowe: zdarzenie, którego zajście leży całkowicie lub częściowo poza zasięgiem
kontroli.
Definuje się:
iloczyn zdarzeń A i B jako równoczesne wystąpienie zdarzenia A i zdarzenia B; (A*B)
1
sumę (alternatywę) zdarzeń A, B, jako wystąpienie zdarzenia A lub zdarzenia B. (A+B)
zdarzenie przeciwne do A zdarzenie zachodzące wtedy gdy A nie zachodzi (~A)
zdarzenie pewne zachodzi zawsze (np. A+(~A));
zdarzenie niemożliwe nie zachodzi nigdy (np. A*(~A)); oznaczamy go symbolem Ć
zdarzenia rozłączne A, B takie, że A*B jest zdarzeniem niemożliwym
(a) (b) (c)
A A A i B B
A lub B
A*B
(A+B)
~
A
nieA
A+(~A)=E
(d) (e)
A B E
Zdarzenie pewne E
(suma wszystkich zdarzeń
możliwych)
Zdarzenia rozłączne A, B
Rys.2. Graficzna ilustracja zdarzeń elementarnych i złożonych:
koła- elementarne zdarzenia losowe, całe ramki wszystkie zdarzenia możliwe
Prawdopodobieństwo zdarzenia liczba wyrażająca stopień możliwości zachodzenia
zdarzenia. Prawdopodobieństwo zdarzenia A czyli P(A) jest równe stosunkowi liczby
przypadków sprzyjających zdarzeniu A (nA) do wszystkich przypadków możliwych (n):
nA
P(A) =
n
Wartość tak zdefiniowanego prawdopodobieństwa ilustrują stosunki pól figur (kół)
reprezentujących zdarzenia A, B na rysunkach powyżej, do pola całej ramki E.
Właściwości prawdopodobieństwa:
1. Jeśli A, B, .. są zdarzeniami rozłącznymi (wykluczają się wzajemnie) to
P(A lub B lub ..)=P(A)+P(B)+.. (patrz rysunek d)
2. Jeśli E jest zdarzeniem pewnym to
P(E)=1 (patrz rysunek e)
Stąd wynika, że dla dowolnego zdarzenia A
0Ł P(A) Ł 1
P(nieA)=1-P(A) (patrz rysunek a)
Dla dowolnych zdarzeń A i B
P(A lub B)=P(A)+P(B)-P(A i B) (patrz rysunki b, c)
Prawdopodobieństwo warunkowe i prawdopodobieństwo całkowite:
2
Mamy dwa zdarzenia losowe A i B. Niech P(B)>0. Jeśli zdarzenia A i B mogą występować
równocześnie to można mówić o prawdopodobieństwie zajścia zdarzenia A pod warunkiem,
że zaszło zdarzenie B, co oznacza się symbolem P(A|B). W tym przypadku zbiór zdarzeń
możliwych redukuje się do zdarzenia B, zatem (zgodnie z rys.2b) mamy:
nAiB nAiB
1 P(A i B) nAiB nB
P(A B) = = = gdzie P(A i B) = ; P(B) =
nB
nB n P(B) n n
n
Prawdopodobieństwo warunkowe nazywane jest prawdopodobieństwem a posteriori (po
uzyskaniu dodatkowej informacji) i na ogół różni się od P(A) zwanego prawdopodobieństwem
a priori (określonym dla dowolnych warunków przy których zachodzi A).
Jeśli w wyniku pewnego doświadczenia losowego realizuje się zawsze jedno z wzajemnie
wykluczających się zdarzeń B1, B2, .. BN (tzn. B1+ B2 +.. BN=E oraz B1*B2=Ć, Bi*Bk=Ć dla
każdej pary zdarzeń Bi, Bk iąk) to dla dowolnego zdarzenia A zachodzi równość:
N
P(A) = P(A B1) P(B1) + P(A B2 ) P(B2 ) + ...P(A BN ) P(BN ) =
P(A i Bn )
n=1
Jest to wzór na prawdopodobieństwo całkowite.
Stosuje się go, gdy prawdopodobieństwa warunkowe P(A|Bn) oraz prawdopodobieństwa P(A)
albo P(Bn) są łatwe do oszacowania lub znane. Prawdopodobieństwo wystąpienia zdarzenia Bn
gdy zaszło zdarzenie A liczy się ze wzoru Bayesa:
P(A Bn ) P(Bn ) P(A Bn ) P(Bn )
P(AiBn )
P(Bn A) = = = ;
N
P(A) P(A)
P(A Bi ) P(Bi )
i=1
Jeśli zdarzenie B nie wpływa na prawdopodobieństwo zdarzenia A to zdarzenia A, B są
zdarzeniami niezależnymi. Wówczas obowiązuje zależność:
P(A B) = P(A); zatem P(A i B) = P(A) P(B)
Zdarzenia A, B są zatem niezależne, gdy mogą występować w różnych okolicznościach, a ich
łączne wystąpienie jest tylko całkowicie przypadkowe. Założenie niezależności zdarzeń jest
często wykorzystywane w obliczeniach probabilistycznych. W ekonometrii mamy na ogół
do czynienia ze zdarzeniami współzależnymi, ale założenie niezależności pozwala dokonać
zgrubnych oszacowań prawdopodobieństw iloczynu zdarzeń.
Zmienne losowe: liczby charakteryzujące rezultat zjawiska losowego
Zmienne losowe dyskretne liczby losowe ze skończonego lub przeliczalnego zbioru
wartości. Na ogół są to liczby całkowite symbolizujące rozważane zdarzenia losowe,
zliczające ich krotność itp.
Zmienne losowe ciągłe: liczby rzeczywiste o losowej wartości, charakteryzujące ilościowo
zjawiska losowe:
Zdarzenia losowe odniesione do liczb losowych dotyczą wystąpienia określonych wartości
zmiennych dyskretnych oraz wystąpienia wartości zmiennych ciągłych w określonych
przedziałach.
Prawdopodobieństwa takich zdarzeń charakteryzują rozkłady prawdopodobieństwa
zmiennych losowych:
3
Dystrybuantą rozkładu zmiennej losowej x nazywamy prawdopodobieństwo wystąpienia
wartości x mniejszej niż argument dystrybuanty (założona wartość zmiennej) X.
F(X) = P(x < X)
Dystrybuanta posiada następujące cechy:
1. F(-Ą)=0;
2. F(Ą)=1;
3. jest funkcją lewostronnie ciągłą i niemalejącą, tzn., jeśli X1
Dystrybuanta zmiennej losowej dyskretnej zmienia się skokowo w punktach
odpowiadających kolejnym wartościom zmiennej.
Rozkład prawdopodobieństwa takich zmiennych wygodniej jest charakteryzować podając
wprost prawdopodobieństwa wystąpienia poszczególnych wartości p(xi). Nazywa się to
krótko rozkładem prawdopodobieństwa zmiennych dyskretnych:
f(x)={p(xi); i=1,2, ...,N}, gdzie N oznacza liczbę możliwych wartości zmiennej x
W przypadku zmiennych losowych ciągłych (dokładnie absolutnie ciągłych patrz
Poradnik [2]) rozkład opisuje się tzw. funkcją gęstości prawdopodobieństwa f(x), którą
definiuje się jako pochodną dystrybuanty względem zmiennej x, tzn. w następujący sposób:
Dx Dx
F(x + ) - F(x - )
d
2 2
f (x) = lim = F(x)
Dxo
Dx dx
Zgodnie z własnością (3) dystrybuanty funkcja f(x) jest nieujemna
Uwaga !! Funkcja gęstości prawdopodobieństwa nie jest prawdopodobieństwem, ale pozwala
obliczyć prawdopodobieństwo wystąpienia wartości X w zadanym przedziale x1, x2 z wzoru:
x2
P(x1 Ł X Ł x2 ) = f (x) dx
x1
X
Wynika stąd, że F(X ) = f (x) dx
-Ą
Ą
oraz f (x) dx = 1
-Ą
Parametry rozkładu prawdopodobieństwa jednowymiarowych zmiennych losowych:
Rozkład prawdopodobieństwa zmiennej ciągłej charakteryzuje się przy pomocy parametrów
zwanych momentami. Moment i-tego rzędu mi(x) definiuje się następująco:
Ą
mi (x) = xi f (x) dx
-Ą
Moment rzędu zerowego jest zawsze równy 1.
Moment rzędu pierwszego zmiennej X nazywa się wartością oczekiwaną zmiennej losowej
X, a operację jego obliczania oznacza się symbolem E(X).
Wartość oczekiwana jest też nazywana wartością przeciętną zmiennej losowej lub nadzieją
matematyczną.
Dla zmiennej losowej ciągłej wartość oczekiwaną wyraża wzór:
4
Ą
def
m1(X ) = E(X ) = X f (X ) dX
-Ą
Dla zmiennej dyskretnej przyjmującej wartości xi z prawdopodobieństwem pi wartość
oczekiwaną oblicza się ze wzoru:
Ą
E(x) = pi xi
i=1
Właściwości wartości oczekiwanej:
1. Każda ograniczona zmienna losowa ma wartość oczekiwaną.
2. Wartość oczekiwana kombinacji liniowej zmiennych losowych jest kombinacją liniową
ich wartości oczekiwanych
E(a1x1 + a2x2 + ... + aN xN ) = a1 E(x1) + a2 E(x2 ) + ... + aN E(xN )
3. Jeśli x1 i x2 są niezależnymi zmiennymi losowymi to
E(x1 x2 ) = E(x1) E(x2 )
Zmienna losowa ciągła x będąca odchyłką zmiennej losowej oryginalnej X od jej wartości
oczekiwanej m1(X) nazywa się zmienną losową scentrowaną:
def
x = X - m1(X );
Oczywiście E(x)=0
Momenty wyższego rzędu można obliczać dla oryginalnych zmiennych lub scentrowanych.
Momenty dla zmiennych scentrowanych nazywa się momentami centralnymi.
Centralny moment rzędu drugiego zmiennej nazywa się wariancją zmiennej s2
x
Ą Ą
s2 = m2[X - E(X )] = m2 (x) = E(x2 ) = - E(X )]2 f (X )dX = x2 f (x)dx
X
[X
-Ą -Ą
Właściwości wariancji:
1. Znając pierwszy i drugi moment oryginalnej zmiennej losowej X można obliczyć jej
wariancję:
2 2
s2 = E(X ) - E2 (X ) = m2 (X ) - m1 (X )
X
bo:
Ą Ą Ą Ą
2 2
X f (X )dX - 2E(X ) X f (X )dX + E (X ) f (X )dX =
[X - E(X )]2 f (X )dX =
-Ą -Ą -Ą -Ą
= m2 (X ) - 2 E2 (X ) + E2 (X ) = m2 (X ) - E2 (X )
2. Jeśli x1, x2 ... xN są niezależnymi zmiennymi losowymi to
2 2 2
E[(a1x1 + a2x2 +K+ aN xN )2 ] = a1 s21 + a2s22 +K+ aNs2
x x xN
Odchyleniem średnim (standardowym) lub dyspersją s zmiennej losowej nazywamy
pierwiastek arytmetyczny z jej wariancji.
Odchyleniem przeciętnym bX zmiennej losowej X nazywamy wartość oczekiwaną modułu
scentrowanej zmiennej x
bX=E(|X-E(X)|)
Parametry pozycyjne rozkładu kwantyle
5
Kwantylem rzędu p zmiennej losowej x nazywamy taką wartość lp zmiennej, że
P(x Ł lp) ł p P(x ł lp) ł 1-p
Kwantyl rzędu p=1/2 nazywa się medianą (jest to wartość zmiennej losowej rozdzielająca jej
zakres na dwie części o jednakowym prawdopodobieństwie wystąpienia p=0.5.
Kwantyle rzędu p=1/4 i nazywają się odpowiednio kwantylem górnym i dolnym.
Kwantyle rzędu p=0.1, 0.2 ....0.9 nazywa się decylami.
Zmienną losową o wartości oczekiwanej 0 i wariancji 1 nazywamy zmienną losową
standaryzowaną.
Jeśli mamy zmienną losową X o wartości oczekiwanej E(X) i dyspersji sX to odpowiadającą
jej zmienną standaryzowaną x uzyskuje się przez przekształcenie
X - E(X )
x =
sX
i odwrotnie, mając zmienną standaryzowaną x, np. odczytaną z tablic rozkładu, uzyskuje się
zmienną X o zadanej : X = x sX + E(X )
Zmienne losowe wielowymiarowe
Jeśli rozważamy kilka zbiorów liczb losowych to mówimy o zmiennej losowej wielowymiarowej.
Dystrybuantę zmiennej wielowymiarowej X=[X1, X2, ... XN ] definiuje się jako
prawdopodobieństwo zdarzenia polegającego na równoczesnym wystąpieniu wszystkich
rozważanych liczb losowych mniejszych od zadanych argumentów dystrybuanty [X1, X2, . XN]
F(X1, X2, ... XN)=P[(x1Wielowymiarowa zmienna losowa ma rozkład absolutnie ciągły jeśli istnieje taka funkcja
f(x1, x2, ... xN ) zwana wielowymiarową gęstością prawdopodobieństwa, że
x1 x2 xN
F(X1, X ,KX ) = f (x1,.x2 ,KxN )dx1.dx2 KdxN
2 N
L
-Ą-Ą -Ą
Zmienne losowe X1, X2, ... XN są niezależne, jeśli ich łączna dystrybuanta jest iloczynem
dystrybuant poszczególnych zmiennych:
F(X1, X2, ... XN)=F(x1Zmienne losowe absolutnie ciągłe są niezależne, jeśli ich wielowymiarowa funkcja gęstości
prawdopodobieństwa jest iloczynem funkcji gęstości dla poszczególnych zmiennych:
Rozkłady brzegowe i warunkowe zmiennych losowych wielowymiarowych
Niech f(x, y) oznacza dwuwymiarowy rozkład zmiennych x, y.
Rozkładami brzegowymi są funkcje:
Ą Ą
fx (x) = f (x, y)dy. f (y) = f (x, y)dx
y
-Ą -Ą
Ą Ą
Oczywiście, obie spełniają warunek podstawowy: fx (x)dx = 1 i f (y)dy = 1
y
-Ą -Ą
6
Rozkładami warunkowymi są natomiast funkcje:
f (x, y0 ) f (x, y0 )
fx (xY = y0 ) = =
Ą
f (y0 )
y
f (x, y0 )dx
-Ą
oraz
f (x0 , y) f (x0, y)
f (y X = x0 ) = =
y
Ą
fx (x0 )
f (x0 , y)dy
-Ą
Kowariancja zmiennych losowych X, Y wartość oczekiwana iloczynu scentrowanych
zmiennych x, y:
Ą Ą
cov(X,Y)=E(xy)=
[X - E(X )][Y - E(Y)] f (X ,Y) dX dY = E(X Y) - E(X ) E(Y)
-Ą-Ą
Współczynnik korelacji rXY zmiennych X i Y to ich kowariancja przeliczona do zakr. [ 1, 1]
rXY=cov(X,Y)/(sxsy)
Współczynnik korelacji przyjmuje wartość 0 gdy zmienne X, Y są niezależne i wartość ą1
gdy są one zależne liniowo (ale tylko liniowo np. X=a+bY, a, b stałe);
Wynika to z następujących rachunków:
cov(XY)=E(aY+bY2)-E(Y)E(X)=aE(Y)+bE(Y2)-E(Y)[a+bE(Y)]=
=b{E(Y2)-[E(Y)]2}=bs2Y
s2X=E[(a+bY)2]-[E(a+bY)]2=E(a2+2abY+b2Y2)-[a+bE(Y)]2=
=a2+2abE(Y)+b2E(Y2)-a2-2abE(Y)-b2[E(Y)]2=b2{E(Y2)-[E(Y)]2}=b2sY2
Zatem sXsY=|b|sY, czyli
rXY=cov(XY)/(sXsY)=b/|b|= ą1
UWAGA !!
Niezerowe, a nawet wysokie wartości współczynnika korelacji dwóch zmiennych
losowych nie oznaczają związku przyczynowo-skutkowego między nimi, a jedynie
współzależność stochastyczną, czyli istnienie wspólnych przyczyn dla obu zjawisk
Procesy stochastyczne
Procesem stochastycznym nazywa się zmienną losową sparametryzowaną czasem (ogólnie
dowolną zmienną skalarną).
Z=(X,t)=Xt
Oznacza to, że każdej chwili czasu to przypisuje się zbiór zmiennych losowych Xto (zwany
zbiorem możliwych realizacji procesu Z w chwili to), z jego wartością oczekiwaną E(Xto) i
rozkładem prawdopodobieństwa f(Xto).
Proces stochastyczny jest zatem szczególnym przypadkiem wielowymiarowej zmiennej
losowej i można dla niego definiować wielowymiarowe rozkłady prawdopodobieństwa
f(Xto, Xt1, .. , XtN ), a także kowariancje cov(Xto, Xt1) odpowiadające różnym to, t1 (zwane
autokowariancjami lub funkcjami korelacyjnymi). Określa się je identycznie jak dla
7
wielowymiarowych zmiennych losowych. Współczynnik korelacji odpowiadający dwóm
różnym wartościom t, nazywa się współczynnikiem autokorelacji procesu.
Ważną klasą procesów stochastycznych są procesy stacjonarne.
Proces stochastyczny jest stacjonarny w węższym sensie jeśli wszystkie jego rozkłady
prawdopodobieństwa nie zależą od czasu, a jedynie od różnic wartości to, t1, .. , tN dla
których są definiowane.
Zatem jednowymiarowe rozkłady prawdopodobieństwa zmiennych Xt dla kolejnych t są
identyczne, czyli f(Xt)=f(X), a rozkłady dwuwymiarowe f(Xt1, Xt2) zależą tylko od różnicy
czasów t=t2-t1, tzn. f(Xt1, Xt2) = f(X, t).
Także autokowariancja i współczynnik autokorelacji zależą tylko od t. Nazywa się je funkcją
korelacyjną KX(t) i funkcją autokorelacji r(t)
def
KX (t) = E(Xt Xt-t ) - (E(X ))2 = E(xt xt-t )
def
E(X X ) - (E(X ))2 E(xt xt-t )
t t-t
rX (t) = =
s2 s2
def
K (t)
X
czyli rX (t) =
s2
gdzie xt oznacza proces scentrowany (tj. E(x)=0).
Wartość funkcji autokorelacji w zerze wynosi zawsze 1, r(0)=1
Proces stochastyczny jest stacjonarny w szerszym sensie jeśli istnieje jego wartość
oczekiwana i jest ona stała w czasie, a funkcje korelacyjne zależą tylko od przesunięcia
czasu t (nie zależą od wartości to, t1):
mX(t)=E(Xt)=mX=const
KX(t1, t2)=E(xt1xt2)=KX(t)
Zatem procesy stochastyczne stacjonarne mają stałą wartość oczekiwaną, a relacje
probabilistyczne między ich wartościami w różnych chwilach czasu są określone
(deterministycznie) przez funkcję autokorelacji r(t).
Mówimy, że funkcja autokorelacji opisuje właściwości dynamiczne procesu
stochastycznego stacjonarnego, natomiast jego właściwości chwilowe (statyczne)
charakteryzuje funkcja gęstości prawdopodobieństwa, czyli rozkład prawdopodobieństwa
Im wolniej maleje funkcja autokorelacji, tym mniej losowe są zmiany w czasie procesu, tzn.
zmiany te są powodowane głównie wewnętrzną inercją procesu, a nie czynnikami losowymi.
Transformata Fouriera funkcji autokorelacji nazywa się funkcją gęstości widmowej mocy
procesu lub spektrum procesu
UWAGA:
Funkcje analogiczne jak funkcja korelacyjna mogą być definiowane dla dwu różnych
procesów stochastycznych stacjonarnych przesuniętych względem siebie w czasie.
Modelowe rozkłady prawdopodobieństwa
Szereg zjawisk losowych można opisać rozkładami prawdopodobieństwa, których gęstości są
stosunkowo prostymi funkcjami analitycznymi zmiennej losowej.
8
Modelowe rozkłady prawdopodobieństwa dla zmiennych dyskretnych (I.E. Brontsztejn,
K.A.Siemeindiajew: Matematyka poradnik encyklopedyczny. Część szósta
Opracowanie danych doświadczalnych,. PWN, Warszawa 1986, str.782)
1. Rozważmy zdarzenie losowe (zwane sukcesem) występujące w pewnym procesie
losowym ze stałym prawdopodobieństwem p>0. Prawdopodobieństwo nie wystąpienia
sukcesu wynosi (1-p).
Jeśli przeprowadzimy n niezależnych doświadczeń, to liczba sukcesów Sn jest liczbą
losową o rozkładzie dwumianowym:
n
ć
P(Sn = k) = pk (1- p)n-k
k
Ł ł
wartość oczekiwana E(Sn)=np; wariancja s2 =np(1-p)
Sn
2. Prawdopodobieństwo łącznej liczby wystąpień X pewnego rzadkiego zdarzenia (o
małym prawdopodobieństwie) oblicza się z rozkładu Poissona, który jest przybliżeniem
rozkładu dwumianowego dla n Ą i p 0 ale tak, że n p l > 0.
lk
P(X = k) = e-l
k!
wartość oczekiwana E(X)= l; wariancja s2 = l
X
Rozkład ten jest w praktyce stosowalny już dla n rzędu kilkudziesięciu, przy l<10.
Przykładem zmiennej X może być liczba klientów zainteresowanych - w pewnym
przedziale czasu - luksusowym artykułem w sklepie odwiedzanym przez wielu klientów
zainteresowanych na ogół innymi artykułami (np. w kiosku). Można go wykorzystać do
oceny opłacalności zamawiania takich artykułów.
3. Rozważmy sytuację jak w (1), ale gdy interesuje nas prawdopodobieństwo zdarzenia
polegającego na wystąpieniu serii k sukcesów, po których następuje brak sukcesu.
Zakładając niezależność kolejnych prób uzyskuje się wyrażenie zwane rozkładem
geometrycznym dyskretnej zmiennej losowej X wyrażającej długość serii niezależnych
sukcesów o jednakowym prawdopodobieństwie wystąpienia:
P(X = k) = pk (1- p)
p p
wartość oczekiwana E(Sn ) = ; wariancja s2 =
X
1- p (1- p)2
Modelowe rozkłady prawdopodobieństwa dla zmiennych ciągłych (I.E. Brontsztejn,
K.A.Siemeindiajew: Matematyka poradnik encyklopedyczny. Część szósta
Opracowanie danych doświadczalnych,. PWN, Warszawa 1986, str.783 i dalsze)
1. Rozkład wykładniczy określa rozstęp czasowy x pomiędzy wystąpieniami zdarzenia
którego prawdopodobieństwo zależy tylko od przedziału czasu w którym się go oczekuje,
a nie zależy od czasu trwania procesu losowego (np. czas pomiędzy nadejściem dwu
kolejnych klientów, czas pomiędzy awariami urządzenia)
0 dla x < 0
a e-ax gdy x ł 0
f (x) = F(x) =
1- e-ax dla x ł 0
0 gdy x < 0
1
wartość oczekiwana E(x)=1/a; wariancja s2 =
X
a2
9
Jest to ciągły odpowiednik rozkładu geometrycznego, gdy x=kdt, a=(1-p)/dt, 1-p=adt,
dt@0. (1-p) prawdopodobieństwo wystąpienia zdarzenia
2. Rozkład równomierny ma zmienna losowa x, gdy może przyjmować z tym samym
prawdopodobieństwem dowolną wartość z przedziału [a-b, a+b], b>0 i nie występuje poza
tym przedziałem:
1 Dystrybuanta i gęstość rozkładu równomiernego
gdy x [a - b, a + b]
f (x) =
2b
0 gdy x [a - b, a + b] 1
F(x)
f(x)
0
x
a+b
a-b
x=a
b2
wartość oczekiwana E(x)=a; wariancja s2 =
X
3
Rozkład równomierny przypisuje się zmiennym losowym, których wartości są naturalnie
ograniczone i wynikają z oddziaływania pewnego czynnika o czysto losowym charakterze
i ograniczonej sile . Przykładowo, taki rozkład może mieć kwota wydawana przez
jednego klienta w małym sklepie spożywczym lub kiosku.
3. Rozkład normalny czyli rozkład Gaussa: dla zmiennej X o wartości oczekiwanej m=E(X)
i dyspersji s (oznaczany symbolem N(m, s)):
ć
1 [X - E(X )]2
f (X ) = exp-
s 2p 2 s2
Ł ł
Dla zmiennej standaryzowanej x rozkład Gaussa N(0, 1) ma postać
ć
1 x2
f (x) = exp-
2
2p
Ł ł
i w tej postaci jest on dostępny w tablicach i generatorach liczb
Dystrybuanta rozkładu Gaussa jest funkcją nieanalityczną zapisywaną w postaci:
x x
2 2
ć ć ć
1 t 1 1 t 1 x
ć
F(x) =
exp- 2 = 2 + 2 2p exp- 2 = 2 erf 2
dt dt 1+
2p
Ł ł
Ł ł Ł ł Ł ł
-Ą -x
gdzie erf(x) (error function funkcja błędu) jest definiowana jako
x
1
2
erf (x) = P(t 2 [-x, x]) = (- )dt
exp t
2p
-x
Wartości funkcji erf(x) dla x=1/ 2 , 2/ 2 , 3/ 2 (czyli w otoczeniu wartości
oczekiwanej o szerokości s, 2s, 3s), wynoszą:
erf(1/ 2 )=0.6827=68.3%
erf(2/ 2 )=0.9545=95.5%
erf(3/ 2 )=0.9973=99.7%
erf(4/ 2 )=0.9999=99.99%
10
Jak widać, zmienne losowe o rozkładzie normalnym praktycznie mieszczą się w zakresie (m-
3s, m+3s) (z prawdopodobieństwem 99.7%). Jest to tzw. zasada trzech sigm.
UWAGA !!! Wartość erf(x) wylicza w MATLABie funkcja o nazwie erf()
f(x)
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
x
0
-3 -2 -1 0 1 2 3
Rys. 3. Rozkład Gaussa N(0,1) dla standaryzowanej zmiennej losowej x
UZUPEANIENIE: Zasada trzech sigm stosuje się dla liczb losowych o dowolnym
rozkładzie.
Ujmuje to ogólnie nierówność Czebyszewa. Niech x oznacza zmienną losową o dowolnym
rozkładzie z ograniczoną wartością oczekiwaną E(x) i ograniczoną wariancją s2. Wówczas
1
P(x - E(x) > k s)Ł
2
k
Zatem, dowolna zmienna losowa mieści się w zakresie trzech sigm z prawdopodobieństwem
co najmniej 90%.
Rozkład normalny mają zmienne losowe, których wartości są zależne od wielu czynników,
przy czym każdy z nich indywidualnie ma mały wpływ na tę wartość. Z taką sytuacją mamy
bardzo często do czynienia w praktyce, w tym również w ekonometrii, dlatego rozkład
normalny odgrywa bardzo ważną rolę w zastosowaniach rachunku prawdopodobieństwa i w
statystyce matematycznej.
Formalnie ujmują to tzw. twierdzenia graniczne (patrz I.E. Brontsztejn, K.A.Siemeindiajew:
Matematyka poradnik encyklopedyczny. Część szósta Opracowanie danych
doświadczalnych,. PWN, Warszawa 1986, Poradnik inżyniera Matematyka Rozdział
XXXII).
11
1. Twierdzenie Lindberga-Leviego: Niech x1, x2, x3, ... xn, będzie ciągiem niezależnych
zmiennych losowych, o jednakowym rozkładzie, posiadających wartość oczekiwaną m i
wariancję s2. Wtedy dla każdego rzeczywistego X spełniona jest relacja:
n
ć
- n m
xk X
ć
1 y2
lim P k =1 < X =
exp- 2 dy
nĄ
s n 2p
Ł ł
-Ą
Ł ł
Oznacza to, że suma n takich zmiennych przy dużej liczbie składników n ma rozkład zbliżony
do normalnego z wariancją równą ns2 i wartością oczekiwaną równą nm.
Można pokazać, że jeśli liczby x1, x2, x3, ... xn mają rozkład równomierny to już dla n
zbliżonych do 10 uzyskuje się praktycznie rozkład normalny.
2. Twierdzenie Lapunowa mówi, że suma zmiennych x1, x2, x3, ... xn zmierza do rozkładu
normalnego także wówczas, gdy mają one różne rozkłady, różne wartości oczekiwane mk i
wariancje s2 , ale muszą mieć odpowiednio silnie ograniczone momenty centralne rzędu
k
trzeciego bk=E(|xk mk|3), tak aby spełniony był następujący warunek:
n
3
bk
k =1
lim = 0 ;
nĄ
s
n
2
gdzie s = jest dyspersją sumy zmiennych.
sk
k =1
Wartość oczekiwana m rozkładu sumy takich liczb jest oczywiście sumą mk dla k=1, 2, ..n.
Z powyższego wynika, że suma niezależnych zmiennych losowych o dowolnych
symetrycznych rozkładach (czyli o zerowych bk, ale posiadających ograniczone wartości
oczekiwane i ograniczone wariancje), zmierza do rozkładu normalnego. Oczywiście,
zbieżność jest w tym przypadku wolniejsza niż dla liczb z Twierdzenia Lindberga-Leviego.
3.Lokalne twierdzenie graniczne Moivre a-Laplace a: Rozkład dwumianowy o ustalonym
prawdopodobieństwie sukcesu p zmierza do rozkładu normalnego, gdy liczba prób n
rośnie do nieskończoności.
Liczbę sukcesów Sn w takiej sytuacji można obliczać jako zmienną losową (ciągłą) o
rozkładzie normalnym z wartością oczekiwaną m=np i wariancją s2= np(1-p). Formalnie
zapisuje się to w postaci:
b
ć
ć
Sn - n p 1 y2
lim Pa Ł Ł b =
exp- 2 dy
nĄ
n p (1- p) 2p
Ł ł
a
Ł ł
Twierdzenie to stosuje się już dla n rzędu kilkudziesięciu. Pozwala ono oszacować przedziały
np. liczby sztuk Sn pewnego artykułu (o stałym popycie) sprzedanych w ciągu dnia, z
wykorzystaniem tablic funkcji błędu erf(x) pakietu MATLAB,
d
gdzie x = ; d - promień przedziału Sn wokół wartości oczekiwanej m=np.
2 n p (1- p)
Przykładowo, dla p=0.02 i n=50 (jeden dzień) mamy E(Sn)=1, sSn=0.99, co daje
P(Sn[0,2])@0.68.
12
Dla n=300 (tydzień) mamy E(Sn)=6, sSn=2.43, co daje P(Sn[4,8])@0.59 oraz
P(Sn[0,12])@0.986, czyli, że sprzedaż tygodniowa praktycznie nie przekroczy 12 sztuk.
W analizie danych wykorzystuje się prawa wielkich liczb.
PRAWO WIELKICH LICZB Chinczyna (Poradnik inżyniera Matematyka str.1072)
Niech x1, x2, x3, ... xn oznacza ciąg niezależnych liczb losowych o jednakowym rozkładzie i
ograniczonej wartości oczekiwanej E(x)=m. Wtedy dla dowolnej liczby dodatniej e zachodzi
równość:
n
ć
1
lim P - m ł e = 0
xk
nĄ
n
k =1
Ł ł
Jeśli dodatkowo niezależne liczby losowe x mają skończoną wariancję to słuszne jest
MOCNE PRAWO WIELKICH LICZB, które mówi, że średnia arytmetyczna ciągu niezależnych
liczb losowych o ograniczonej wartości oczekiwanej i wariancji zmierza z
prawdopodobieństwem 1 do ich wartości oczekiwanej.
n
ć 1
Plim = E(x) = 1
xk
ŁnĄ n k =1 ł
4. ELEMENTY STATYSTYKI MATEMATYCZNEJ
Literatura:
1. J.Greń: Statystyka matematyczna modele i zadania, PWN, Warszawa 1982
2. "Ekonometria. Metody i analiza problemów ekonomicznych". Pod red. K. Jajugi; Wydawnictwo
AE we Wrocławiu, Wrocław 1999.
3. I.E. Brontsztejn, K.A.Siemeindiajew: Matematyka poradnik encyklopedyczny. Część
szósta Opracowanie danych doświadczalnych,. PWN, Warszawa 1986
4. Poradnik inżyniera Matematyka Rozdział XXXIII
5. Edward Nowak: Zarys metod ekonometrii zbiór zadań, PWN, Warszawa 1994
6. John Freund Podstawy nowoczesnej statystyki, PWE, Warszawa 1968
Pojęcia podstawowe statystyki (wg J.Greń: Statystyka matematyczna)
Metody statystyki stosuje się w ekonometrii w celu badania właściwości probabilistycznych
zmiennych ekonometrycznych. Zmienne te reprezentują pewne cechy badanej zbiorowości,
czyli populacji generalnej. Poddawane analizie dane (liczby losowe) są widziane jako losowa
próba populacji generalnej. Aby wynik badania był miarodajny, należy zadbać o
reprezentatywność próby. Uzyskuje się ją przez odpowiedni dobór (losowanie, rejestrację)
elementów próby z populacji generalnej.
Losowe wartości cechy w próbie o liczności n traktuje się jako jedną wartość n-wymiarowego
wektora losowego. Zbiór wszystkich możliwych wartości tej cechy w próbie o liczności n
nazywa się przestrzenią próby.
Rozkład populacji, to rozkład wartości badanej cechy w populacji generalnej. Zwykle zakłada
się, że rozkład populacji jest zbliżony do pewnego rozkładu modelowego (patrz poprzedni
rozdział).
13
Na podstawie zebranych danych losowych można obliczać statystyki z próby, czyli dowolne
funkcje zebranych zmiennych losowych. Statystyki są także zmiennymi losowymi.
Takimi statystykami mogą być zależności pozwalające na oszacowanie parametrów
rozkładów prawdopodobieństwa badanych cech populacji generalnej. Nazywa się je
estymatorami parametrów rozkładów. Wynikiem zastosowania estymatora jest estymata
poszukiwanego parametru rozkładu.
Estymator nieobciążony estymator a pewnego parametru a spełniający równość E(a)= a,
co oznacza, że estymator szacuje wartość a bez błędu systematycznego, a więc pozwala
znalezć faktyczną wartość parametru.
Estymator zgodny estymator a pewnego parametru a spełniający warunek:
lim P( an - a < e) = 1
nĄ
tzn, estymator, który jest stochastycznie zbieżny do wartości parametru. Gdy używa się
estymatora zgodnego, to zwiększanie liczności próby zmniejsza błąd estymacji.
Ogólnie, estymator zgodny może być obciążony, a estymator nieobciążony może nie być
zgodny.
Estymator efektywny estymator o możliwie małej wariancji.
Badania statystyczne zmiennych ekonometrycznych obejmują:
1. analizę losowości zmiennych i określenie ich rozkładu prawdopodobieństwa;
2. obliczanie czyli estymację parametrów tych rozkładów;
3. weryfikację hipotez statystycznych dotyczących rozkładów populacji generalnej
Analiza losowości zmiennych i określenie ich rozkładu prawdopodobieństwa
Podstawowe estymatory rozkładów zmiennych losowych
1. Zgodnie z prawem wielkich liczb, jeśli cecha x populacji generalnej jest zmienną losową o
ograniczonej wariancji s2 i ograniczonej wartości oczekiwanej m, a ciąg wartości x1, x2, ...xn
tej cechy w próbie jest ciągiem niezależnych liczb losowych, to estymatorem zgodnym,
nieobciążonym i najefektywniejszym wartości m=E(x) jest jej średnia arytetyczna:
n
-
1
E(x)=m @ x =
xk
n
k =1
Średnia arytmetyczna liczb o rozkładzie normalnym ma rozkład normalny o wartości
oczekiwanej m (takiej jak liczby uśredniane) i wariancji
s
s2 =
xs
n
Jeśli uśredniane liczby mają rozkład inny niż normalny to rozkład średniej arytmetycznej
zmierza (dla n zmierzającego do nieskończoności) do rozkładu normalnego N(m,s n)
(patrz Twierdzenie Lindberga-Leviego)
1. Zgodnie z prawem wielkich liczb, jeśli cecha x populacji generalnej jest zmienną losową o
ograniczonej wariancji s2 i ograniczonej wartości oczekiwanej m, a ciąg wartości x1, x2, ...xn
tej cechy w próbie jest ciągiem niezależnych liczb losowych, to estymatorem zgodnym,
nieobciążonym i najefektywniejszym wartości m=E(x) jest jej średnia arytetyczna:
n
-
1
E(x)=m @ x =
xk
n
k =1
Średnia arytmetyczna liczb o rozkładzie normalnym ma rozkład normalny o wartości
oczekiwanej m (takiej jak liczby uśredniane) i wariancji
14
s
s2 =
xs
n
2. Estymatorem zgodnym, nieobciążonym i najefektywniejszym wariancji cechy x na
podstawie ciągu wartości x1, x2, ...xn tej cechy w próbie o własnościach jak wyżej, jest
średniokwadratowa odchyłka od wartości średniej obliczana wg wzoru:
n n
- -
1 1 n
2
E[(x - m)2 ] = s2 @ s2 = - x)2 = - (x)2
(xk xk
n -1 n -1 n -1
k =1 k =1
Występująca w mianowniku różnica (n-1) powoduje, że estymator jest nieobciążony. Wynika
ona z faktu, że dane w liczbie n zostały już wykorzystane do obliczenia wartości średniej.
Zatem liczba jeszcze nie wykorzystanych danych wynosi (n-1).
Formalnie można to wykazać wychodząc z zależności:
n n
- -
1 1 n
2
E[(x - m)2 ] = E{ - x)2} = } - E{(x)2}
(xk E{xk
M M M
k=1 k=1
gdzie M jest nieznanym jeszcze dzielnikiem estymatora.
Ze wzorów podanych w poprzednim rozdziale wynika, że:
2
s
2 2 2
x
E{xk } = s + m2 , E{x2} = s + m2 = + m2
x xs
n
Po podstawieniu do wzoru wyżej mamy:
2 2
n
1 s ć
n (n -1) s (n -1)
2 2 2
x x
s = E[(x - m)2] = s + m2 - - m2 = = s
ć x x
M n M n m
k=1
Ł ł Ł ł
Wynika stąd, że aby uzyskać tożsamość, M powinno być równe (n-1).
Dzielenie przez n daje zatem błąd estymatora (obciążenie) jakkolwiek pozostaje on zgodny.
Dla dużych n można ten błąd zaniedbać i dzielić przez n.
n n
- -
1 1
2
s2 @ - x)2 = - (x)2
(xk xk
n n
k =1 k =1
Estymata wariancji obliczona jak wyżej dla liczb o rozkładzie normalnym ma rozkład c2
(chi-kwadrat) patrz I.E. Brontsztejn, K.A.Siemeindiajew: Matematyka poradnik
encyklopedyczny. Część szósta.
(
Zmienna losowa h ma rozkład c2 o k stopniach swobody, gdy gęstość phn) (c) wyraża
się wzorem
k
k
-
-1
2
2
e c
ć
(
phk ) (c) =
k
2
k Ł ł
2
2 Gć
2
Ł ł
Suma kwadratów N niezależnych liczb losowych o rozkładzie normalnym N(0,1) ma rozkład
c2 o N-1 stopniach swobody. Zatem
(N -1)s2 2 (N -1)s2 b ( -1)
P < s < = phN (c)dc = a
ż
b a
a
4.Estymatorem kowariancji dwóch zmiennych losowych x, y jest wyrażenie:
n
- -
1
E[(x - mx ) (y - my )] = cov(x, y) @ Kxy = - x) (yk - y)
(xk
n -1
k =1
lub w przybliżeniu (dla dużych n)
15
n n
- - - -
1 1
Kxy @ - x) (yk - y) = yk - x y
(xk xk
n n
k =1 k =1
Zachodzi równość: Kxy=Kyx
Analogicznie, dla wielowymiarowych zmiennych losowych X=[X1, X2, .. XM], gdzie X1, X2, ..
XM są wektorami kolumnowymi (o jednakowej długości) losowych wartości kolejnych cech w
próbie, liczy się macierz kowariancji K wg wzoru macierzowego:
T
- -
1
T
K = X X - X X
n
-
gdzie T oznacza transpozycję macierzy, X jest wektorem wierszowym wartości średnich
- -
T
kolejnych wektorów X1, X2, .. XM, , zatem X (wektor kolumnowy) pomnożony przez X
(wektor wierszowy) daje macierz kwadratową o wymiarze MxM, czyli taką jak macierz K.
Jeśli wektory X1, X2, .. XM są wektorami wierszowymi (takie są na ogół produkowane przez
procedury pakietu MATLAB), to macierz kowariancji K liczy się wg wzoru:
T
- -
1
T
K = X X - X X
n
Macierz kowariancji K jest symetryczna, tzn KT=K
5.Estymator funkcji autokowariancji i autokorelacji
Jeśli ciąg liczb losowych jest uporządkowany wg czasu ich rejestracji, z jednakowym
rozstępem czasowym Dt między próbkami, czyli do kolejnych próbek można przypisać
czas liczony numerem kolejnym próbki, to ciąg jest reprezentacją dyskretną procesu
stochastycznego (albo inaczej zdyskretyzowanym procesem stochastycznym).
Jeśli badany proces stochastyczny jest stacjonarny, to estymatorem jego funkcji
autokowariancji z rozstępem m (czyli z rozstępem czasowym t=mDt) ma postać
n
- -
1
Rxn (m) = - x) (xk -m - x)
(xk
n -1- m
k=m+1
Funkcję autokowariancji liczy się na ogół dla wystarczająco długich ciągów, aby można było
zastosować wzór:
2
n n
- - -
1 1
ć
Rxn (m) @ - x) (xk -m - x) = xk xk - x
(xk -m
n - m n - m
Ł ł
k =m+1 k =m+1
Jeśli ciąg x1, x2, ...xn jest ciągiem niezależnych liczb losowych reprezentujących stacjonarny
proces stochastyczny, to
s2 dla m = 0
x
lim Rxn (m) =
n-Ą
0 dla m > 0
Estymatorem unormowanej funkcji autokowariancji (czyli autokorelacji) rxn(m) jest
wyrażenie
Rxn (m)
rxn (m) =
2
sx
16
Dyskusje ocen autokowariancji ciągu rnm dla kolejnych m zawiera podręcznik
[Mańczak Nahorski] str.66. Dla dużych n i dużych m mamy:
Ą
1
cov(rm, rm+1) = rk+1
rk
n
k=-Ą
rm prawdziwa wartość funkcji
Stąd wynika oszacowanie wariancji rm
Ą
1 ć
2
var(rm ) = 1 + 2 (i)
rk
n
Ł k=1 ł
Jeśli proces xi jest liniowy i spełnia warunki:
Ą Ą Ą Ą
2 2
xn = xn-k < Ą ak |< Ą
ak ak | | k | ak < Ą
k=-Ą k=-Ą k=-Ą k=-Ą
to wielowymiarowy rozkład zmiennych n(rm - rm ) N(0,W ) .
nĄ
Dla dużych n istotność funkcji autokorelacji można badać metodą Boxa:
1. Zakładamy, że r1, ..., rk są zerowe
2. obliczamy wariancje rm ze wzoru (i) i sprawdzamy hipotezę j.w. dla m=1
3. jeśli należy ją odrzucić przyjmujemy r1=r1 i sprawdzamy hipotezę dla m=2
4. w taki sam sposób sprawdzamy dla dalszych m
Szeregi rozdzielcze i estymacja rozkładu prawdopodobieństwa
Histogram
a) dzielimy przedział zmienności zbioru liczb x na pewną liczbę m podprzedziałów o
szerokości Dxi. i grupujemy liczby w klasy takie, aby w i-tej klasie znalazły się wszystkie
- -
liczby xk ( x - 0.5Dxi , x + 0.5Dxi ] Należy przy tym tak dobrać liczbę klas lub ich
i i
szerokości, aby w każdej klasie znalazło się co najmniej 10 liczb Uzyskuje się w ten
sposób szereg rozdzielczy zmiennej x.
b) dla każdego przedziału liczymy wysokość słupka histogramu pi wg wzoru:
ni 1
pi =
n Dxi
- -
Wartość pi jest estymatą prawdopodobieństwa Pi{x ( x - 0.5Dxi , x + 0.5Dxi ] }
i i
-
c) robimy wykres słupkowy wartości pi względem x i uzyskujemy histogram.
i
Na tle histogramu wskazane jest wykreślenie przypuszczalnego rozkładu teoretycznego (lub
zestawu rozkładów) populacji generalnej, z parametrami E(x) oraz s2 obliczonymi jak wyżej.
Jeśli danych jest wystarczająco dużo, najwygodniej jest przyjąć stałą szerokość klas szeregu
rozdzielczego. Można do tego wykorzystać funkcję hist(x,k) z pakietu MATLAB. Tworzy
ona szereg rozdzielczy wektora x z liczbą przedziałów k, a następnie wykreśla histogram dla
liczności ni próbek w podprzedziałach (a nie dla pi jak zalecono w punkcie b)
Histogram można wykonać dla dowolnego zbioru liczb, ale ma on sens, gdy liczby te są
liczbami losowymi. Sam kształt histogramu nie rozstrzyga kwestii losowości danych,
jakkolwiek może sugerować, że są lub nie są one losowe. Hipotezę losowości należy przyjąć
wcześniej na podstawie analizy przyczyn losowości danych, obecności pewnych oddziaływań
deterministycznych, czy wreszcie przypuszczalnego rozkładu prawdopodobieństwa danych.
Na ogół dla zmiennych ciągłych rozkład ten winien być zbliżony do normalnego, w pewnych
17
przypadkach do równomiernego. Stwierdzenie dużych rozbieżności histogramu od takiej
hipotezy może być spowodowane:
a) niejednorodnością próbki, tzn. występowaniem kilku (np. dwóch) klas danych o różnych
rozkładach;
b) zależnością parametrów rozkładu (głównie wartości oczekiwanej) od pewnych czynników
nielosowych (np. czasu lub zmiennych egzogenicznych)
W przypadku (a) zbiór danych należy rozdzielić na podzbiory jednorodne statystycznie.
W przypadku (b) należy znalezć najpierw odpowiednie zależności ekonometryczne badanej
zmiennej (endogenicznej) od znanych czynników egzogenicznych, a następnie poddać
analizie probabilistycznej odchyłki losowe danych od tych zależności. Będzie to przedmiotem
następnego rozdziału opracowania.
Estymacja przedziałowa
Jak, wspomniano, estymaty są liczbami losowymi o określonym rozkładzie
scharakteryzowanym dyspersją i wartością oczekiwaną. W związku z tym wynik estymacji
można przedstawić jako estymatę wartości oczekiwanej parametru oraz przedział, w którym
z zadanym prawdopodobieństwem zwanym poziomem ufności mieści się prawdziwa
wartość parametru. Przedział taki nazywa się przedziałem ufności, sam sposób postępowania
estymacją przedziałową.
Zwykle podaje się przedział ufności o szerokości 1s, 2s lub 3s. Dla rozkładu normalnego
estymatora odpowiada to poziomowi ufności odpowiednio: 68.3%, 95.5% i 99.7%
Zgodność histogramu z założonym rozkładem prawdopodobieństwa testuje się m.in.
przy pomocy statystyki c2 (kryterium c2 Pearsona).
Niech F(x) oznacza założoną (teoretyczną) dystrybuantę zmiennej losowej x, f(x) funkcję
gęstości prawdopodobieństwa, D1, D1, ..., Dj, ..., DL ciąg rozłącznych przedziałów histogramu
o wartości średniej xj, nj - liczbę danych w j-tym przedziale. Wówczas, dla NĄ statystyka
2
ć
n
j
dF(x) - N 2
L L
(N f (x )D - n )
Dj
j j j
2 ł
c = N @
Ł
N f (x )D
j=1 j=1
j j
dF(x)
Dj
ma rozkład c2 o L-1 stopniach swobody. W praktyce wystarcza, aby minj N f(x)Dj>10.
5. Modele ekonometryczne
Literatura:
1. Manikowski A., Tarapata Z.: Prognozowanie i symulacja rozwoju przedsiębiorstwa.
WSE Warszawa 2002
2. Pawłowski Z. " Zasady predykcji ekonometrycznej" PWN, Warszawa 1982.
3. Zeliaś A. "Teoria prognozy" PWE, Warszawa 1997.
4. Dittmann P., Metody prognozowania sprzedaży w przedsiębiorstwie, wyd. 6, Wyd. AE
Wrocław, 2002.
5. Gajda J.B., Prognozowanie i symulacja a decyzje gospodarcze, C.H.Beck Warszawa,
2001.
6. Radzikowska B. (red.), Metody prognozowania. Zbiór zadań, wyd. 3, Wyd. AE Wrocław,
2001.
18
7. Zeigler B.P, Teoria modelowania i symulacji, PWN Warszawa, 1984.
8. K. Molenda, M. Molenda, Analiza i prognozowanie szeregów czasowych, Placet,
Warszawa 1999
9. E. Nowak. (red.) Prognozowanie gospodarcze. Metody, modele, zastosowania,
przykłady. Placet 1998
10. Cieślak M. Prognozowanie gospodarcze. Wydawnictwo AE Wrocław, 1998.
11. Henry Theil: Zasady ekonometrii, PWN, Warszawa, 1979
12. Zbigniew Pawłowski: Ekonometria, PWN, Warszawa 1969
6. G.E.P.Box, G.M.Jenkins: Analiza szeregów czasowych, PWN, Warszawa, 1983
Modelem ekonometrycznym nazywa się zależność stochastyczną wartości oczekiwanych
zmiennych endogenicznych od deterministycznych zmiennych egzogenicznych.
Schemat ujęcia zagadnienia zilustrowano na Rys.1.
Model ekonometryczny zapisuje się ogólnie w postaci:
^
y = f (X ) ;
gdzie X jest wektorem wartości zmiennych, przyjętych jako istotnie oddziałujące na zmienną
objaśnianą y. Zmienna objaśniająca y ma wartość:
^
y = y+ e
gdzie e jest błędem modelu spowodowanym z założenia tylko czynnikami losowymi.
Jeśli wektor X zawiera tylko zmienne egzogeniczne (zewnętrzne) przypisane do tej samej
chwili czasu co zmienna endogeniczna, to model jest równaniem algebraicznym opisującym
proces ekonomiczny w stanie wewnętrznej równowagi, a model nazywa się zależnością
statyczną.
Modele statyczne buduje się w celu znalezienia ilościowych i stałych w czasie powiązań
zmiennych egzogenicznych z endogenicznymi.
Jeśli co najmniej jedna spośród zmiennych X jest zmienną endogeniczną zarejestrowaną w
poprzednich chwilach czasu, lub zmienne egzogeniczne brane są z wcześniejszych chwil
czasu niż y, to model jest dyskretnym zapisem równania różniczkowego (czyli równaniem
różnicowy) i opisuje przebieg zjawiska dynamicznego (uwarunkowanego wewnętrzną
dynamiką badanego procesu ekonomicznego). Taki model nazywa się modelem dynamiki
zjawiska. Modele takie wykorzystuje się w ekonometrii głównie do prognozowania przebiegu
zjawiska w przyszłości (matematyczne prognozowanie procesów ekonomicznych, np. cen,
zapotzrzebowania itp.)
W ekonometrii czynniki losowe odgrywają zwykle dużą rolę, stąd konieczność poszukiwania
raczej prostych zależności f(X).
W praktyce najczęściej stosuje się modele jednoczynnikowe, gdzie mamy tylko jedną
zmienną objaśniającą. Opracowanie wiarygodnych modeli wieloczynnikowych jest zadaniem
bardzo trudnym i wymaga często długotrwałych badań.
Szczególnym przypadkiem modeli jednoczynnikowych są takie, gdzie jedyną zmienną
objaśniającą jest czas t. Zależność ekonometryczną tego typu nazywa się często trendem, a
wyliczanie odchyłek danych od takiej zależności ekstrakcją trendu. Formalnie, zależność
taka jest zależnością statyczną, niemniej może być (i często jest) wykorzystywana do analizy
19
dynamicznego przebiegu zjawiska. W szczególności stosuje się takie modele do
prognozowania zjawisk (tzw. modele Browna).
Modele regresyjne, liniowe
Bardzo ważną grupę modeli ekonometrycznych stanowią modele liniowe o ogólnej postaci:
K
^
yn = ukn (X ) = Un A
ak
k =0
gdzie Un jest wektorem wierszowym tzw. wejść uogólnionych [u0n, u1n, ..., uKn], A
wektorem kolumnowym nieznanych współczynników modelu A=[a0, a1, ..., aK].
Wejścia uogólnione są to funkcje algebraiczne zmiennych objaśniających, zadane wraz ze
wszystkimi koniecznymi współczynnikami, tak aby mając wartości zmiennych Xn w
kolejnych próbkach n=1..N, można było obliczyć wszystkie wartości Un dla n=1, ..N.
Zatem, dane o wejściach X pozwalają obliczyć całą macierz U złożoną z wierszy Un.
^ ^
Niech Y, Y oznaczają wektory kolumnowe wartości yn , yn dla n=1, .. N; E wektor
kolumnowy kolejnych wartości zakłóceń e1, e2, ..., eN. Można ogólnie zapisać formuły:
^ ^
Y = U A + E ; lub Y = Y + E gdzie Y = U A
W ekonometrii najczęściej wykorzystuje się takie modele bądz jako modele wieloczynnikowe
liniowe, tj. biorąc Un Xn, lub jako modele jednoczynnikowe, w których wejścia u(x) są
prostymi funkcjami, przeważnie jednomianami. W ostatnim przypadku mamy modele
wielomianowe (najpopularniejsze).
Postać wejść U oraz ich liczba to struktura modelu regresyjnego. Struktura modelu
regresyjnego musi być ustalona arbitralnie. Mając tę strukturę należy obliczyć parametry A.
Wyznaczanie parametrów A nazywa się identyfikacją modelu. Realizuje się ją w następujący
sposób:
a) Dokonuje się pomiarów w liczbie N>K+1, tak aby nadmiar danych w stosunku do liczby
nieznanych współczynników pozwolił dobrać optymalne współczynniki
b) Definiuje się kryterium jakości modelu
c) Oblicza się wartości współczynników rozwiązując zadanie optymalizacji polegające na
minimalizacji kryterium (b)
Najprostszy algorytm uzyskuje się przyjmując jako kryterium (b) sumę kwadratów błędów
modelu. Nazywa się to metodą najmniejszych kwadratów MNK.
Takie zadanie ma rozwiązanie analityczne, tzn. optymalny wektor współczynników modelu
wyraża się wzorem macierzowym:
-1
T T
A = [U U] U Y
UWAGA !!!
Wyznaczone w powyższy sposób parametry A mają sens tylko wówczas, gdy macierz U jest
dobrze uwarunkowana, tzn. wynik jej odwracania jest słabo zależny od błędów
numerycznych. Sprawdza się to procedurą svd() z biblioteki MATLAB. Uwarunkowanie nie
zależy od wartości wyjść procesu (i ich składnika losowego) ale tylko od struktury funkcji
regresji i zmienności wejść (im mniejsza zmienność tym gorsze uwarunkowanie).
Uwarunkowanie można zatem poprawić tylko przez zmianę struktury modelu lub zebranie
bardziej zróżnicowanych danych wejściowych.
20
Warto podkreślić, że taki sam wzór określa współczynniki optymalnej aproksymacji dowolnej
funkcji y(u) zadanej w N punktach. Model ekonometryczny ma być jednak zależnością
stochastyczną, a nie aproksymatą pewnej funkcji deterministycznej. Oznacza to, że błedy
modelu chcemy interpretować jako spowodowane tylko czynnikami losowymi, a nie
wynikające z arbitralnego doboru wartości funkcji.
Aby uzyskany model mógł być interpretowany jako zależność stochastyczna winny być
spełnione następujące warunki (podane przez Gaussa) (patrz Z.Pawłowski: Ekonometria):
a) zmienne objaśniające X, a więc także U winny być nielosowe (czyli dokładnie znane), i U
nie mogą być liniowo współzależne (ogólnie macierz U musi być dobrze
uwarunkowana, co jest konieczne, aby obliczenia dawały sensowne rezultaty)
b) składnik losowy e zmiennej objaśnianej musi mieć zerową wartość oczekiwaną E(e)=0, i
skończoną oraz stałą wariancję (niezależną od czasu przypisanego do kolejnych danych)
c) ciąg {en, n=1,2, ..N}musi być ciągiem niezależnych liczb losowych
d) składnik losowy e nie może być skorelowany ze zmiennymi objaśniającymi
uwzględnionymi w modelu.
Jeśli te założenia są spełnione to estymatory A współczynników modelu są zgodne i
nieobciążone, a ich macierz kowariancji wyraża się prostym wzorem:
T 2
KA = [U U]-1 sr
2
gdzie sr oznacza estymatę wariancji reszt modelu.
Pozwala to wyznaczać nie tylko samą funkcję regresji, ale także przedziały ufności dla tej
2
funkcji oraz prognoz zmiennej y, regresji w oparciu o oceny wariancji funkcji regresji sym i
2
wariancji przewidywanej zmiennej objaśnianej sy , które oblicza się ze wzorów:
2 T 2 2 2
sym = un KA un ; sy = sym + sr
Model jednoczynnikowe wielomianowe mają wejścia uogólnione o postaci:
u0n 1; (stała modelu)
u1n= xn
u2n= (xn)2
.........
uKn= (xn)K
Model z ograniczeniami równościowymi [Pawłowski str.120]
BA=C
B znana macierz, C znany wektor, A wektor współczynników oryginalnych.
Zmodyfikowany wektor parametrów wyraża się wzorem:
A* = [I - HB]A + H C
gdzie
H = [UT U]-1 BT [B(UT U)-1BT ]-1
Macierz kowariancji ma postać:
T
2
KA = [I - HB][UT U]-1[I - HB] sr
*
Funkcje korelacyjne a zależności regresyjne
W przypadku szeregów czasowych, stosowaną powszechnie analizę współczynnika
korelacji zmiennych losowych można łatwo rozszerzyć, badając korelacje szeregów
21
przesuniętych względem siebie w czasie, czyli funkcje korelacyjne. Dla szeregów XNd={x1, ..
xN-d}, YN={yd+1, .. yN}przesuniętych o d próbek definiuje je formuła:
def def N
E{(YN - ysr)(XN -d - xsr)} 1 1
(1)
Ryxd = @ R* = - ysr)(xn-d - xsr)
yxd (yn
N - d sysx n=d +1
E{(YN - y)2}E{(XN -d - x)2}
N N
1 1
, (2)
sx = - xsr)2 sy = - ysr)2
(xn-d (yn
N - d N - d
n=d +1 n=d +1
gdzie sx, sy, , , ysr, xsr oznaczają dyspersje, wartości oczekiwane i średnie ciągów XNd i YN.
y x
Wzór (1) można zastosować dla dowolnych ciągów, ale funkcja Ryxd jest miarodajna, gdy
są one stacjonarnymi szeregami czasowymi [6].
Przyjmijmy, że szeregi YN, XNd mają zerowe wartości średnie (ysr=0, xsr=0) (co można
zawsze uzyskać prowadząc analizy dla szeregów scentrowanych) i są powiązane liniową
funkcją regresji:
, (3)
yn = axn-d + zn wn = axn-d
gdzie zn są próbkami zakłócenia, - wartością oczekiwaną E{yn} (zależną liniowo od xn-d).
wn
Jeśli zakłócenia zn mają rozkład normalny z zerową wartością oczekiwaną oraz stałą wariancją
(niezależną od n), i nie są skorelowane z poprzednimi wartościami (E{znzn-m}=0 dla m >0) i ze
zmienną x (E{znxn-d}=0), to optymalny, zgodny i nieobciążony estymator współczynnika a
wyraża się wzorem wynikającym z metody najmniejszych kwadratów [Paw]:
N
sssy sy
1 (4)
a = xn-d = R* = R*
yn
N 2
sx yxd sx yxd
2
n=d +1
xn-d
n=d +1
Dyspersję se błędów zależności (3) wyraża wzór:
def N N
sy
1 1
(5)
se = - axn-d )2 = (yn - R* xn-d )2 = sy 1- (R* )2
(yn yxd
N - d N - d sx yxd
n=d +1 n=d +1
a dyspersja współczynnika a obliczonego ze wzoru (4) wynosi:
sy 1- (R* )2
1 se 1 (6)
yxd
sa = se = =
N
sx N - d sx N - d
2
xn-d
n=d +1
Miarą zasadności stosowania modelu regresyjnego (3) jest statystyka Studenta (t) jego
współczynnika a:
def
a N - d
(7)
ta = = R*
sa yxd 1- (R* )2
yxd
Wynika stąd też zależność między statystyką Studenta i współczynnikiem korelacji, którą
wykorzystuje się do testowania istotności korelacji [A.Manikowski, Z.Tarapata:
Prognozowanie i symulacja rozwoju przedsiębiorstwa. WSE Warszawa 2002 (str. 172)]:
ta 2
(7a)
R* =
yxd
N - d - 2 + ta 2
Dla odpowiednio długich, nieskorelowanych ciągów (N-d >30) statystyka ta ma rozkład
Gaussa N(0,1) [Pa]. Uzyskanie wartości ta>3 oznacza zatem, że zależność (3) jest istotna na
poziomie istotności 0.3%. Zatem, jeśli Ryxd ą 0, to wykazanie tego wymaga ciągu o długości
(8)
N - d > 9((R* )-2 -1)
yxd
22
Bardziej praktyczną miarą przydatności modelu jest jednak stopień redukcji niepewności
informacji o wartościach yn uzyskiwanych wg wzoru (3) względem oceny trywialnej, jaką
daje średnia wartość ysr szeregu YN . Zgodnie ze wzorem (5) miarę tę wyraża formuła:
se
(9)
= 1- (R* )2
yxd
sy
Jeśli szeregi YN, XN-d są stacjonarne, to dla odpowiednio dużych wartości N, np.
spełniających relację (8), uzyskuje się @Ryxd =const, sy@sy, sx@sx, se@se oraz:
R*
yxd
2
, (10)
se /s = 1- Ryxd
y
gdzie sx, sy, se są stałymi odchyleniami standardowymi zmiennych x, y i błędów modelu (3),
Wówczas właściwości prognoz uzyskiwanych przez zastosowanie modelu (3), t.j.
(11)
wn+d = axn
można miarodajnie ocenić wg wzoru (9), a istotność zależności (3) wg. wzoru (8).
W przeciwnym wypadku ekstrapolacja zawarta w formule predyktora (11) może wprowadzać
znaczący błąd, większy niż wynikałoby to ze wzoru (9), szczególnie dla małych N.
Modele regresyjne dynamiki procesów (modele ARMAX)
Ważną rolę odgrywają regresyjne, liniowe modele dynamiki procesów wynikające z równań
różnicowych opisujących przebieg procesu w dyskretnych chwilach czasu:
I J K
yn = yn-i + xn- j-m + zn-k + zn
ai b j g j
i=1 j=1 k =1
gdzie pierwszą sumę nazywa się członem autoregresyjnym (AR), a sumy druga i trzecia są
średnimi ruchowymi (moving average MA) dla zmiennych egzogenicznych X (np. sterowań
branych z zadanym opóznieniem m próbek) oraz zakłóceń z (zakłocenia z są z założenia
ciągiem niezależnych liczb losowych o zerowej wartości oczekiwanej) .
Modele takie omawia obszernie monografia G.E.P.Box, G.M.Jenkins: Analiza szeregów
czasowych, PWN, Warszawa, 1983.
W pakiecie MATLAB są procedury do identyfikacji współczynników modeli ARMAX.
W przypadku, gdy wejścia są nieznane (np. gdy interesują nas przebiegi wejść zakłócających
widzianych jako efekt czynników wyłącznie losowych), stosuje się model ARMA:
I K
yn = yn-i + zn-k + zn
ai g j
i=1 k =1
Jest on często stosowany do prognozowania takich zakłóceń (jako alternatywa dla model
Browna).
UWAGA !!! W modelu ARMA (a także ARMAX) zakłada się, że zakłócenia z mają
wartość oczekiwaną zerową. Jeśli nie jest to spełnione, a w szczególności, gdy
wartość ta jest zmienna, lepiej zastosować model ARMA dla różnic wyjść y i
zakłóceń. Taki model nazywa się modelem ARIMA (scałkowany model ARMA)
Z drugiej strony, pomijając drugą sumę modelu ARMAX uzyskuje się model dynamiki
procesu dla przypadku, gdy mamy słabe zakłócenia zewnętrzne, a z reprezentuje szumy
pomiarowe. Taki model identyfikuje się metodą analizy regresji, biorąc jako wejścia
uogólnione u wektor: un=[yn-1, yn-2, ... yn-I, xn-1-m, xn-2-m, ... xn-J-m ]
oraz jako wejście ciąg: yn=[yn, yn-1, ... yn-I+1].
23
Wektor współczynników: A=[a1, a2, ...aI, b1, b2, ...bJ];
Procedura identyfikacji może dalej przebiegać jak dla modeli statycznych, jakkolwiek lepsze
rezultaty daje metoda zmodyfikowana, zwana metodą zmiennych instrumentalnych (IV). Jest
ona omówiona np. w podręczniku: A.Niederliński: Komputerowe systemy sterowania.
W pakiecie MATLAB jest procedura do identyfikacji współczynników modelu tą metodą
(funkcja iv()).
Modelowanie charakterystyk statycznych metodą analizy regresji
(Uzupełnienie [Duda 2003])
Modele statyczne buduje się w celu znalezienia ilościowych i stałych w czasie powiązań
zmiennych wejściowych z wyjściowymi. Jeśli czynniki losowe odgrywają dużą rolę, zaleca
się poszukiwanie raczej prostych zależności f(U).
Bardzo ważną grupę modeli statystycznych procesów stanowią zależności liniowe (ze
względu na nieznane parametry) o ogólnej postaci:
K
^
yn = jkn (Un ) = Fn A (1.6.3)
ak
k =0
gdzie jn jest wektorem wierszowym tzw. wejść uogólnionych [j0n, j1n, ..., jKn], A
wektorem kolumnowym nieznanych współczynników modelu A=[a0, a1, ..., aK].
Wejścia uogólnione są to funkcje algebraiczne zmiennych objaśniających, zadane wraz ze
wszystkimi koniecznymi współczynnikami, tak aby mając wartości zmiennych Un w
kolejnych próbkach n=1..N, można było obliczyć wszystkie wartości jn dla n=1, ..N. Zatem,
dane o wejściach U pozwalają obliczyć całą macierz F złożoną z wierszy jn. Zakres
zmienności poszczególnych kolumn macierzy F wyznacza N-wymiarowe pole korelacji
modelu.
^ ^
Niech Y, Y oznaczają wektory kolumnowe wartości yn , yn dla n=1, .. N; Z wektor
kolumnowy kolejnych wartości zakłóceń z1, z2, ..., zN:
^ ^
Y = F A + Z ; lub Y = Y + Z gdzie Y = F A (1.6.4)
W sterowaniu najczęściej wykorzystuje się takie modele bądz jako modele wieloczynnikowe
liniowe, tj. biorąc jn=f(Un), lub jako modele jednoczynnikowe, w których wejścia jn(u) są
prostymi funkcjami (przeważnie jednomianami) jednej zmiennej wejściowej.
Postać wejść F oraz ich liczba to struktura modelu regresyjnego, która musi być
ustalona arbitralnie. Mając tę strukturę należy obliczyć parametry A. Wyznaczanie
parametrów A nazywa się identyfikacją modelu. Jest ona realizowana w następujący sposób:
d) rejestruje się pomiary wejść i wyjść w liczbie N>K+1, tak aby nadmiar danych w stosunku
do liczby nieznanych współczynników pozwolił dobrać optymalne współczynniki A
e) definiuje się kryterium jakości modelu
f) oblicza się wartości współczynników rozwiązując zadanie optymalizacji polegające na
minimalizacji kryterium (b),
g) przy pomocy testów statystycznych usuwa się z modelu mało istotne składniki.
Najprostszy algorytm identyfikacji uzyskuje się przyjmując jako kryterium (b) sumę
kwadratów błędów modelu. Nazywa się to metodą najmniejszych kwadratów MNK. Takie
zadanie ma rozwiązanie analityczne, tzn. optymalny wektor współczynników modelu wyraża
się wzorem macierzowym ([Mań79], [MaN81], [Gór97]):
-1
A = [FT F] FTY (1.6.5)
24
Wzór ten można zapisać w postaci splotowej:
A = GY (1.6.5a)
def
-1
gdzie G =[FTF] FT jest macierzą, której wiersze są odpowiedziami impulsowymi filtru
przekształcającego ciąg Y na współczynniki jego aproksymaty minimalnokwadratowej.
Zostanie wykazane w dalszej części tego rozdziału, że suma elementów wierszy k=0, 1, ....K
macierzy G spełnia równości:
N N
= 1 oraz = 0 dla k=1, ...K (1.6.5b)
gon gkn
n=1 n=1
Należy zwrócić uwagę, że wyznaczone wg wzoru (1.6.5) parametry A mają sens tylko
wówczas, gdy macierz F jest dobrze uwarunkowana, tzn. wynik jej odwracania jest słabo
zależny od błędów numerycznych (patrz rozdz.3.2). Uwarunkowanie nie zależy od wartości
wyjść procesu (i zakłóceń) ale tylko od struktury funkcji regresji i zmienności wejść (im
mniejsza zmienność tym gorsze uwarunkowanie). Można je zatem poprawić tylko przez
przeskalowanie (normalizację) wejść uogólnionych lub zmianę struktury modelu, lub zebranie
bardziej zróżnicowanych danych wejściowych. Dokładność obliczeń można również
zwiększyć stosując tzw. faktoryzację pierwiastkową i przekształcenie Hausholdera [Nied85].
Wzór (1.6.5) jest ogólną formułą minimalno-kwadratowego dopasowania danych do
założonego modelu, a więc nie zawsze daje wyniki będące estymatami w sensie
statystycznym. Określa bowiem również współczynniki optymalnej aproksymacji dowolnej
funkcji y(u) zadanej w N punktach. Dla niektórych zastosowań poszukiwana zależność (1.6.3)
może spełniać jedynie rolę aproksymaty pewnej złożonej funkcji w polu korelacji.
Przykładem są zależności wykorzystywane do iteracyjnej optymalizacji on-line omawiane w
rozdz.3.5.3.2. W wielu przypadkach wymaga się jednak, aby model był adekwatny również w
pewnym otoczeniu pola korelacji. Funkcja (1.6.4) musi być wówczas zależnością
stochastyczną co oznacza, że błędy modelu winny być powodowane tylko czynnikami
losowymi, bez znaczącego wpływu arbitralnego doboru struktury modelu.
Model (1.6.3) o współczynnikach obliczonych wg wzoru (1.6.5) jest zależnością
stochastyczną, jeśli spełnione są następujące warunki podane przez Gaussa (patrz [Paw69]):
a) zmienne wejściowe U, a więc także F muszą być nielosowe (czyli dokładnie
znane), a macierz F musi być dobrze uwarunkowana,
b) składnik losowy z zmiennej wyjściowej musi mieć zerową wartość oczekiwaną
E{z}=0, i skończoną oraz stałą wariancję (niezależną od czasu przypisanego
do kolejnych danych)
c) ciąg {zn, n=1,2, ..N} musi być ciągiem niezależnych liczb losowych
d) składnik losowy z nie może być skorelowany ze zmiennymi wejściowymi
uwzględnionymi w modelu, tzn. E{FT Z} = 0 .
Jeśli spełnione jest założenie (a), to błąd losowy współczynników wyraża się wzorem:
dA = [FT F]-1FT Z (1.6.6)
Wynika stąd, że przy spełnionym założeniu (d) estymator (1.6.5) jest nieobciążony, tzn.:
E{dA}= 0, a macierz kowariancji współczynników modelu wyraża ogólnie wzór:
def
T
KA = E{dA(dA)T } = [FT F]-1FT E{Z Z }F[FT F]-1 = [FT F]-1FT KZF[FT F]-1 (1.6.7)
gdzie KZ oznacza macierz autokowariancji zakłóceń, która przy założeniach (b) i (c) wynosi
2 2 2
Is . Estymatę zgodną sr wariancji s zakłóceń (wariancji reszt modelu) wyraża wzór:
z z
25
2
N
^
1
2
sr =
ć yn - yn (1.6.8)
N - K -1
Ł ł
n=1
Zatem, jeśli spełnione są założenia (a-d), to estymator (1.6.5) współczynników modelu jest
zgodny i nieobciążony, a macierz kowariancji współczynników wyraża się prostym wzorem:
2
KA = [FT F]-1 sr (1.6.9)
Jeśli dodatkowo zakłócenia losowe mają rozkład Gaussa, to wzór (1.6.5) daje
najefektywniejsze estymaty współczynników modelu o przyjętej strukturze [Mań79].
Wzór (1.6.8) daje estymator zgodny wariancji resztowej. Można to wykazać, podobnie
jak dla estymatora wariancji, wychodząc z zależności:
N N
1 1
2 2 2
E{sr } = E{ yn - wn )2} = (E{yn} - E{wn}) (1.6.9a)
(
M M
n=1 n=1
gdzie M jest nieznanym jeszcze dzielnikiem estymatora.
Ze wzorów podanych w rozdziale o rachunku prawdopodobieństwa wynika, że:
2 2
~n
E{yn} = sz + y2 (1.6.9b)
2 2
~n
E{wn} = s + y2 (1.6.9c)
yn
2
~n
gdzie y = E{yn}. s oznacza wariancję modelu.
yn
Pierwszy składnik po prawej stronie wzoru (1.6.9c) jest wartością oczekiwaną estymaty
wariancji modelu. Błąd modelu dla n-tej obserwacji wejść uogólnionych jn wynosi (zgodnie z
1.6.5b i 1.6.6):
N K N
dwn = jn[FTF]-1FT Z = gki zi = zi (1.6.9d)
jnk cni
i=1 k=0 i=1
2
Estymatę s2(Un ) wariancji modelu s dla obserwacji Un można obliczyć na dwa sposoby.
y yn
Pierwszy, to wykorzystanie macierzy kowariancji współczynników (1.6.9):
T
s2(Un ) = jn (Un )KAjn (Un ) (1.6.9e)
y
2
Jeśli KZ= Is , to jej wartość oczekiwana wynosi:
z
2 T 2 2
E{s2 } = s = {jn[FTF]-1jn }s = g sz (1.6.9f)
yn yn z n
def
T
gdzie g = jn[FTF]-1jn są współczynnikami rozkładu błędu generowanego przez
n
zakłócenia dla poszczególnych obserwacji jn będących wierszami macierzy F.
Przykładowe wartości współczynników filtrujących gkn i rozkładu gn dla 100 obserwacji
modelu wielomianowego 3. stopnia pokazano na rys.1.
Wykażemy, że
N
= K +1 (1.6.9g)
g n
n=1
Niech fij oznacza ij-ty element macierzy F=[FTF]-1. Współczynnik gn wyraża wzór:
26
f00 f01 f0K
ł ł ł
ę ś ę ś ę ś
10 11 1K
g = [jn0,jn1,..,jnK]ę:f śjn0 + [jn0,jn1,..,jnK]ę:f śjn1 + ... + [jn0,jn1,..,jnK]ę:f śjnK =
n
ę ś ę ś ę ś
ę ś ę ę ś
fK fK1ś fKK
0
f00 f01 f0K
ł ł ł
ę ś ę ś ę ś
10 11 1K
= [jn0jn0,jn1jn0,....,jnKjn0]ę:f ś + [jn0jn1,..,jnKjn1]ę:f ś + ... + [jn0jnK ,.,jnKjnK]ę:f ś
ę ś ę ś ę ś
ę ś ę ę ś
fK fK1ś fKK
0
co można zapisać w skrócie:
f0k
ł
ę ś
K
1k
g = jnk,...,jn1jnk,..,jnKjnk]ę:f ś (1.6.9g0)
n [j
n0
ę ś
k=0
ę
fKkś
Zsumowanie elementów gn daje:
f0k
ł
N K N N N
łę f1k ś
= jnk, ..., jnk,.., jnk śę: ś (1.6.9g1)
g n jn0 jn1 jnK
ę
n=1 k =0 n=1 n=1 n=1 ę ś
ę ś
fKk
Zwróćmy uwagę, że wektor wierszowy jest tu k-tym wierszem macierzy FTF, natomiast
wektor kolumnowy to k-ta kolumna macierzy odwrotnej F, tj. [FTF]-1. Iloczyn skalarny
takich wektorów musi być zawsze równy 1, co wynika z oczywistej relacji [FTF][FTF]-1=I.
Zatem powyższa suma wynosi zawsze K+1, co jest wykazaniem zależności (1.6.9g) !!!!.
Wynika stąd ogólne twierdzenie algebry macierzowej:
Ślad macierzy U*[UTU]-1UT jest zasze równy rzędowi macierzy U
W taki sam sposób wykazuje się zależność (1.6.5b).
Sposób drugi, to zastosowanie wzoru (1.6.9d).
N K N
dwn = jn[FTF]-1FT Z = gki zi = zi
jnk cni
i=1 k=0 i=1
Po podstawieniu wzorów (1.6.9b) i (1.6.9cb) do wzoru (1.6.9a) i uwzględnieniu (1.6.9f) oraz
(1.6.9g) mamy:
2
N N N
1 s ć (N - K -1)
2 2 2 2 2 2
z
~n yn ~n 1 ć z
E{sr } = (s + y2 -s - y2)= Ns - = s
= N -
z s yn g n z
M M M M
n=1 Ł n=1 ł Ł n=1 ł
Wynika stąd, że aby uzyskać estymator nieobciążony, M powinno być równe (N-K-1).
Wzór (1.6.9) pozwala wyznaczyć przedziały ufności dla funkcji regresji oraz prognoz
2
zmiennej y, w oparciu o oceny wariancji funkcji regresji sym i wariancji przewidywanej
2
zmiennej objaśnianej sy , które oblicza się ze wzorów:
27
2 2 2
sy = sym + sr (1.6.10)
Odpowiedzi impulsowe MNK dla K=3, W=I Rozklad bledow modelu w polu korelacji gn(n)
0.16
0.14
0.14
0.12
0.12
0.1
0.1
0.08
0.06
0.08
0.04
0.02
0.06
0
0.04
-0.02
-0.04
20 40 60 80 100 20 40 60 80 100
n=1, ..., 100 n=1, ..., 100
Uzyskanie dobrych wyników modelowania, szczególnie poza polem korelacji,
wymaga trafnego doboru funkcji regresji, tj. przekształceń jn(U) generujących wejścia
uogólnione, a także przyjęcia odpowiedniej liczby członów (K+1) modelu (1.6.3). Dla
procesów wielowejściowych dobór przekształceń j jest na ogół trudnym problemem,
wymagającym dużego doświadczenia, intuicji i wiedzy technologicznej. Jeśli rozbieżności
między funkcją regresji pierwszego rodzaju a modelem są duże, to formuła (1.6.9) jest
nieadekwatna, a co za tym idzie oszacowania przedziałów ufności są niemiarodajne (na ogół
zbyt optymistyczne). Z drugiej strony, zastosowanie zbyt złożonej funkcji regresji prowadzi
często do złego uwarunkowania zadania identyfikacji, lub powoduje bardzo niekorzystne
właściwości modelu poza polem korelacji, ujawniające się zwykle szybkim wzrostem
wariancji modelu obliczanej wg (1.6.10). Jest to typowy efekt nadparametryzacji modelu,
której skutkiem jest nadmierna wrażliwość wyniku estymacji na zakłócenia.
Podobne zależności można uzyskać dla uogólnionej MNK.
A = [FTWF]-1FTWY (1.6.7)
def
T T
KA = E{dA(dA)T } = [FTWF]-1FTWE{Z Z }W F[FTWF]-1 =
(1.6.7)
T
[FTWF]-1FTWKZW F[FTWF]-1
Jeśli jako macierz wagową przyjmiemy WT=a(KZ)-1, gdzie a jest dowolną stałą, to
współczynniki A nie zależą od a, natomiast macierz KA wyraża się wzorem:
KA = a[FTWF]-1 (1.6.9h)
2
W przypadku, gdy KZ jest diagonalna macierz W=diag( a /s ) . Obliczymy sumę ważonych
zn
kwadratów błędów modelu:
2
^
2
ć
N N
^
ć
ETWE = yn - yn = a
wn yns- yn (1.6.8h)
Ł ł
n=1 n=1
zn
Ł ł
Podobnie jak we wzorach (1.6.8-1.6.9f) obliczymy wartość oczekiwaną tej sumy:
28
n
k,n
g
g
2
^
2
ć
2 2 2
N N N
ć s
ć
yn
zn
(1.6.9i)
E{ETWE} = a
yns- yn } = a E{yn} - E{wn} = as
2 2 2 2
s zn - zn
s s
n=1 n=1 n=1
zn Ł zn ł
Ł ł
Ł ł
2
Estymatę s2(Un ) wariancji modelu s dla obserwacji Un można obliczyć wg macierzy
y yn
kowariancji współczynników (1.6.9h):
T
s2(Un ) = jn (Un )KAjn (Un ) (1.6.9j)
y
W rozważanym przypadku jej wartość oczekiwana wynosi:
2 T - T 2
E{s2 } = s = a{jn[FTWF]-1jn }= jn[FT KZ1F]-1jn = g szn (1.6.9k)
yn yn n
def
T
gdzie g = jn[FTF]-1jn są współczynnikami rozkładu błędu generowanego przez
n
zakłócenia dla poszczególnych obserwacji jn, takimi samymi jak dla MNK wzór (1.6.9f).
Zatem, zgodnie ze wzorem (1.6.9g) mamy
2
N N
ć
s
yn
= (K + 1) (1.6.9l)
= g n
2
s
n=1 n=1
zn
Ł ł
Zatem
E{ETWE} = a(N - K -1) (1.6.9m)
a więc estymatorem zgodnym stałej a jest
2
N
^
ETWE 1
ć
a = = yn - yn (1.6.9n)
wn
(N - K -1) (N - K -1)
Ł ł
n=1
Macierz kowariancji zakłóceń należy obliczać ze wzoru:
2
N
^
1
ć
KA = [FTWF]-1 yn - yn (1.6.9o)
wn
(N - K -1)
Ł ł
n=1
Przykładowe współczynniki gkn i rozkład błędów dla K=3 i N=100 pokazano na rysunku.
Odpowiedzi impulsowe MNK dla K=3, W=I*1/n Rozklad bledow modelu w polu korelacji gn(n)/wn
11
0.5
10
0.45
9
0.4
8
0.35
7
0.3
6
0.25
5
0.2
4
0.15
3
0.1
2
0.05
1
0
20 40 60 80 100 20 40 60 80 100
n=1, ..., 100 n=1, ..., 100
W celu uzyskania modelu o najkorzystniejszych właściwościach należy
wyspecyfikować w miarę liczny zbiór wejść uogólnionych, a następnie wybrać z niego
możliwie mało liczny podzbiór wejść istotnych w sensie statystycznym. Suboptymalne
procedury wyboru zwane regresją krokową [MaN81] (ang. step-wize regression) są
stosunkowo łatwe do algorytmizacji. Najprostsza i zwykle wystarczająco skuteczna jest
procedura odrzucania [MaN81], oparta na iteracyjnym testowaniu hipotezy Ho o nieistotności
kolejnych współczynników aj modelu:
29
k,n
n
n
g
g
/w
Ho: E{aj}=0; H1: E{aj}ą0;
W typowym przypadku, gdy zakłócenia z mają rozkład normalny, testowanie hipotezy Ho
opiera się na statystyce t Studenta [Mań79], z założonym poziomem istotności testu a.
Funkcja gęstości prawdopodobieństawa rozkładu Studenta i związek poziomu istotności a z
wartością krytyczną tkr statyki mają postać [Gór97]:
l +1
Gć
tkr
1 2 1
Ł ł
f (t,l) = 1- f (t,l)dt = a (1.6.11)
l+1
l
pl
2
-tkr
2
Gć ć
t
2
Ł ł
1+ l
Ł ł
gdzie l>0 jest parametrem całkowitym (liczbą stopni swobody), G() - funkcją gamma Eulera.
Niech saj oznacza estymatę dyspersji współczynnika aj, obliczonego w m-tej iteracji, w
której model (1.6.3) zawierał (Km+1) członów (saj.= K ). Jeśli hipoteza Ho jest prawdziwa,
Ajj
to iloraz tj= aj/saj ma rozkład Studenta t o liczbie stopni swobody l=N-Km-1 [Mań79]. Zatem,
rozpoczynając od nadmiarowego (ale zapewniającego dobre uwarunkowanie zadania)
zestawu członów równania (1.6.3), w kolejnej iteracji oblicza się współczynniki A według
wzoru (1.6.5) oraz wartości statystyki tj= aj/saj dla każdego członu obecnego w modelu. Jeśli
spełniona jest relacja:
| t |= min | ti |< tkr (1.6.12)
j
i
to przyjmuje się hipotezę Ho dla j-tego członu równania, bo nie ma podstaw do jej odrzucenia
(oznacza to założenie aj=0, a więc usunięcie j-tego członu z równania), po czym rozpoczyna
się następna iteracja. Procedura kończy się, gdy relacja (1.6.12) nie jest spełniona, co oznacza,
że wszystkie składniki modelu są istotne na poziomie istotności a.
1.6.3. Modele regresyjne dynamiki procesów
W klasycznym ujęciu punktem wyjścia do identyfikacji dynamiki procesu SISO jest model
dyskretny w postaci ARMAX (CARMA):
I J K
yn = yn-i + un- j-m + zn-k + zn (1.6.13)
ai b j g j
i=1 j=1 k =1
gdzie zakłócenia z są z założenia ciągiem niezależnych liczb losowych o zerowej wartości
oczekiwanej.
Modele takie omawiają obszernie monografie [BoJ83], [MaN81] i [Nied85]. Procedury ich
identyfikacji zawarte są w pakietach wspomagania projektowania układów automatyki (np.
MATLAB).
W przypadku, gdy wejścia są nieznane (np. gdy interesują nas przebiegi wejść zakłócających
widzianych jako efekt czynników wyłącznie losowych), stosuje się model ARMA:
I K
yn = yn-i + zn-k + zn (1.6.14)
ai g j
i=1 k =1
Jest on często stosowany do prognozowania zakłóceń.
Należy zaznaczyć, że w modelu ARMA (a także CARMA) zakłócenia z mają z założenia
wartość oczekiwaną zerową. Jeśli nie jest to spełnione, a w szczególności, gdy wartość ta jest
zmienna, lepiej zastosować model ARMA dla różnic wyjść y i zakłóceń. Taki model nazywa
się modelem ARIMA (scałkowany model ARMA)
Z kolei, pomijając drugą sumę modelu ARMAX uzyskuje się model dynamiki procesu
dla przypadku, gdy mamy słabe zakłócenia zewnętrzne, a z reprezentuje szumy pomiarowe.
Taki model może uwzględniać także opóznienie d. Identyfikuje się je metodą analizy regresji,
biorąc jako wejścia uogólnione j wektor:
30
jn=[yn-1, yn-2, ... yn-I, un-1-m-d, un-2-m-d, ... un-J-m-d ] (1.6.15)
oraz jako wyjście ciąg:
yn=[yn, yn-1, ... yn-I+1]. (1.6.16)
Wektor współczynników ma postać:
A=[a1, a2, ...aI, b1, b2, ...bJ] (1.6.17)
Procedura identyfikacji może dalej przebiegać jak dla modeli statycznych, z macierzą F
zbudowaną z wierszy (1.6.15), ale wzór (1.6.5) daje obciążone estymaty współczynników A,
gdyż zakłócenia zn modelu (1.6.14) są skorelowane ze zmiennymi wejściowymi yn-i.
Lepsze rezultaty można uzyskać przez zastosowanie zmodyfikowanej procedury identyfikacji,
zwanej metodą zmiennych instrumentalnych (IV) [Nied85]. W celu wyznaczenia
nieobciążonej estymaty współczynników A modelu dynamicznego definiuje się macierz W
zmiennych instrumentalnych, taką aby W nie była skorelowana z zakłóceniami zn, a
równocześnie była silnie skorelowana z wejściami uogólnionymi j, tj. :
1 1
T T
lim W Z = 0 , lim W F ą 0 (1.6.18)
N Ą N Ą
N +1 N +1
Symbol lim oznacza tu granice stochastyczne [Nied85], a druga macierz musi istnieć i być
dobrze uwarunkowana. Jeśli znana jest taka macierz W, to nieobciążonym i zgodnym
estymatorem współczynników A jest formuła:
T T
A = [W F]-1W Y (1.6.19)
Macierz W konstruuje się iteracyjnie, wychodząc od macierzy F i zamieniając występujące w
niej wartości zakłócone yn-i ich estymatami uzyskanymi w poprzedniej iteracji. Można
wykazać, że jeśli zakłócenia wyjścia obiektu są szumem białym, to procedura taka jest
zbieżna [Nied85].
Dla obiektów z opóznieniem wartość d dobiera się metodą prób tak, aby zminimalizować
wariancję resztową modelu.
Model dyskretny (1.6.14) winien być przekształcony do postaci ciągłej. Jest to możliwe,
jeśli wszystkie pierwiastki rzeczywiste wielomianu o współczynnikach (1, a1, a2, ...aI) są
dodatnie. Jeśli występują pierwiastki ujemne, to model nie ma odpowiednika ciągłego, co
może być spowodowane zbyt małą częstotliwością próbkowania, pominiętą nieliniowością,
lub zbyt słabym pobudzeniem u obiektu w stosunku do zakłóceń.
Od lat 70. rozwijane są także metody identyfikacji modeli ciągłych ([UnR87],
[FBy84]). Punktem wyjścia jest równanie różniczkowe zwyczajne o postaci:
I J K
y = y(i) + u( j) + z(n+k ) + z (1.6.20)
ai b j g j
i=1 j=0 k =1
które splata się obustronnie z odpowiedzią impulsową g(t) pewnego filtru określonego na
nośniku zwartym [0, T], tj. g(t)ą0 dla t[0, T] oraz g(t)= 0 dla t[0, T].
Jeśli funkcja g(t) posiada ciągłe pochodne do J-tej włącznie dla t[0, T], to wykorzystując
właściwości splotu uzyskuje się model regresyjny w postaci:
I J K
^
Y = ji + j + jI +J +k +1 (1.6.21)
ai b j j+I +1 gk
i=1 j=0 k =1
gdzie
T T
def
^ def
(i)
Y = g(t)y(t - t)dt ; ji = g (t)y(t - t)dt; dla i=1,...., I; (1.6.22)
0 0
T T
def def
( j) (k )
j = g (t)u(t - t)dt dla j=0, .....J; jI +J +k +1 = g (t)z(t - t)dt (1.6.23)
j+I +1
0 0
31
Obliczając numerycznie wartości powyższych splotów dla odpowiednio długiego ciągu
dyskretnych wartości t, uzyskuje się wektor Y i macierz F pozwalające obliczyć
współczynniki modelu ciągłego przy pomocy wzoru (1.6.5). Udoskonalone wersje tej metody
omawiają prace [FBy84], [ByF92].
W prognozowaniu zakłóceń stochastycznych, jako alternatywę dla modeli ARMA
można wykorzystać jednoczynnikowe modele statyczne, w których jedyną zmienną
wejściową jest czas t (tzw. modele Browna). Zależność tego typu nazywa się często trendem,
a wyliczanie odchyłek danych od takiej zależności ekstrakcją trendu [DuK99a].
Szczególne znaczenie dla syntezy klasycznych układów regulacji mają proste modele
niskiego rzędu z opóznieniem. Są one wyznaczane przez aproksymację odpowiedzi obiektu
na wymuszenie skokowe, modelem liniowym o wymaganej strukturze. Identyfikacja takich
modeli może być realizowana technikami regresji statycznej. Jako model procesu przyjmuje
się teoretyczną odpowiedz układu, a zmienną niezależną jest czas. Formalnie, prowadzi to do
nieliniowego zadania optymalizacji. Jednakże w szczególnie interesującym przypadku, gdy
poszukuje się optymalnego modelu pierwszego rzędu z opóznieniem, problem można
sprowadzić do zadania liniowo-kwadratowego poprzez przekształcenie logarytmiczne
wyjścia. Odpowiedz skokowa takiego obiektu ma bowiem postać:
(t - d) ł
w(t) = K expć-
ę(1- T ś (1.6.24)
Ł ł
gdzie K oznacza wzmocnienie, T stałą czasową, d - opóznienie.
Stosując przekształcenie logarytmiczne, wzór (1.6.24) sprowadzamy do postaci:
(t - d)
ln(K - w(t))= ln(K) - (1.6.25)
T
W wyniku eksperymentu uzyskujemy ciąg Y={y1, y2, yp, ......, yN, ....., ye, ....., yf} wartości
wyjść zarejestrowanych w chwilach czasu U = {t1, t2, tp, ..., tN, .., te, .., tf} (niekoniecznie z
jednakowym rozstępem). Wartość te jest czasem po którym wyjście nie wykracza poza
przedział K(1ą e) do końca okresu obserwacji tf, gdzie e jest odpowiednio małą liczbą
dodatnią (np. 0.02). Wartość tf winna być istotnie większa niż te (około 5 krotna wartość
przewidywanej stałej czasowej). Jeśli wyjście jest silnie zakłócone, to jako e można przyjąć
oszacowanie błędu standardowego zakłócenia, w przedziale (te, tf) przekroczenia zakresu
K(1ą e) winny być tylko losowe, a ich liczba rzędu 1/3 liczby danych w tym przedziale.
Uśredniając dane {ye, ....., yf} oblicza się wzmocnienie K, a następnie wektor wyjść GpN
modelu dla chwil tn z zakresu [tp tN]:
def
GpN ={ln(K - y(tn )): n = p,....N} (1.6.26)
Oczywiście, przekształcenie obejmuje również zakłócenia zn, będące różnicą między błędem
estymacji K i zakłóceniami chwilowymi wyjścia. W związku z tym wartość tN musi być
odpowiednio mniejsza niż te. Zgodnie ze wzorem (1.6.25) elementy gn ciągu GpN można
wyrazić zależnością liniową:
gn =a+b un + zn (1.6.27)
def def
d 1
gdzie a = ln(K) + , b = - (1.6.28)
T T
Wartości współczynników a i b uzyskuje się metodą analizy regresji liniowej ze wzoru
(1.6.27), a następnie oblicza się parametry d i T ze wzoru (1.6.28).
Czas te, jak również czas początkowy tp należy dobrać tak, aby błędy liniowej aproksymacji
ciągu GpN nie były duże, np. tak aby zminimalizować wariancję resztową modelu, a
zakłócenia resztowe zn nie wykazywały istotnego spadku ze wzrostem tn.
32
ZESTAWIENIE FUNKCJI PAKIETU MATLAB, UŻYTECZNYCH W EKONOMETRII
abs(x) wartość bezwzględna zmiennej x lub długość wektora x
arx([y u],[na nb nd]) identyfikacja wg.najmniejszych kwadratów modelu ARX ciągu y z
wejsciem x; na rząd mianownika, nb rząd licznika, nd - opóznienie
arxmax([y u],[na nb nc nd]) identyfikacja wg.najmniejszych kwadratów modelu ARMAX
ciągu y z wejsciem x; na rząd mianownika, nb rząd licznika, nc rząd modelu
MA zakłócenia losowego, nd - opóznienie
bar(x,y) rysowanie wykresu słupkowego (histogramu) elementów wektora y na pozycjach
wskazanych przez wartości wektora x.
chi2inv(p, lst_swob) wartość chi2 na poziomie istotnosci p przy lst_swob stopni swobody
cdf() - różne funkcje rozkładu.
cov(x) wariancja wektora x; macierz kowariancji, gdy x jest macierzą, której wiersze są
wektorami zmiennych losowych
erf(x) - funkcja błędu, tj. prawdopodobieństwo wystąpienia zmiennej losowej o rozkładzie
normalnym w przedziale E(x)ąx2
floor(x) obcięcie liczy x do wartości całkowitej
hist(y,n) rysowanie histogramu szeregu rozdzielczego elementów wektora y w n
jednakowych przedziałach; wywołanie: [y,x]=hist(y,n); powoduje tylko policzenie
wektorów histogramu do narysowania funkcją bar(x,y);
inv(A) odwrócenie macierzy A (kwadratowej)
iv([y u],[na nb nd],nf,mf) identyfikacja morelu ARX metodą zmiennych instrumentalnych
mean(x) wartość średnia wektora x
median(x) mediana wektora x
a=polyfit(x,y,n) obliczenie współczynników wielomianu przybliżającego
minimalnokwadratowo zależność wektora y od x; n rząd wielomianu; a(1), a(2),
... współczynniki wielomianu od najwyższej potęgi a(1) do stałej zapisanej jako a(n)
rand generator liczb losowych (niezależnych) o rozkładzie równomiernym z zakresu [0, 1]
randn - generator liczb losowych (niezależnych) o rozkładzie normalnym standaryzowanym,
czyli N(0,1)
std(x) odchylenie standardowe (dyspersja) elementów wektora x
svd(A) dekompozycja singularna macierzy A (sprawdzanie dopuszczalności odwracania
macierzy)
xcorr(x) podaje (m.innymi) funkcje autokorelacji wektora x (nieunormowaną)
Dodatkowe informacje o funkcji uzyskuje się komendą MATLABa:
help nazwa
33
Teskt programu do badania rozkładów zmiennych (Rys.3-6)
% ============= Przetwarzanie danych
ekonometrycznych==============
clear min t tt syg ws hsyg wsyg wsred wsred2 wfilt wfilt2 wsred2sh
wfilt2sh aa okno;
Tf=10; Tsr=2*Tf; Tsrd=Tsr; dt=1;
ldan=256; ldan=ldan*2+Tsrd; Sig=1;
normalny=1;
% -------- generowanie sygnalu o rozkl.rownomiernym ----------
if(normalny==0)
for(i=1:ldan) syg(i)=Sig*(rand-0.5)*sqrt(12); end
else
% ------ generowanie sygnalu o rozkl.normalny ---------------
for(i=1:ldan) syg(i)=Sig*randn; end
end
% =============== Koniec przygotowania danych ===============
% ============== Przetwarzanie wlasciwe
=============================
% ----- Symulacja sygnalow niskoczestotliwosciowych wsred ---
Tsr=Tsrd,; Sum=0; ld=1; sygf=syg(1); clear okno;
for(i=1:ldan)
Sum=Sum+syg(i);
if(i>Tsr) Sum=Sum-okno(ld); ldusr=Tsr; else ldusr=ld; end;
okno(ld)=syg(i);
if(ld==Tsr) ld=0; end;
ld=ld+1;
wsred(i)=Sum/ldusr;
end
Ldan=ldan;
ww=syg(Tsr:Ldan); clear syg; syg=ww; ldan=length(syg); clear ww;
ww=wsred(Tsr:Ldan); clear wsred; wsred=ww; clear ww;
% ============ rysunki ===================
for(niejednor=0:1)
if(niejednor==1)
for(i=1:ldan) wsred(i)=1+10*i/ldan+syg(i); end
syg(ldan/2:ldan)=3*syg(ldan/2:ldan);
end
% ------- Normalizacja i centrowanie sygnalu ---------------
sum=0; Sum=0; for(i=1:ldan) sum=sum+wsred(i); Sum=Sum+syg(i); end;
sred=sum/ldan; Srsyg=Sum/ldan;
for(i=1:ldan) t(i)=i-1; wsred(i)=wsred(i)-sred; syg(i)=syg(i)-
Srsyg; end;
Sig=std(syg); Sigs=std(wsred);
% ----------------------------------------------------------
figure; subplot(2,1,1)
plot(t,syg,'y',t,wsred,'b')
axis([0 max(t) min([min(syg) min(wsred)]) max([max(syg)
max(wsred)])])
subplot(2,1,2)
plot(t,wsred,'b')
axis([0 max(t) min(wsred) max(wsred)]);
xlabel(sprintf('Usrednianie'))
34
% ======= histogram =====================
podpis='abcd'; kk=1;
figure;
for(typ=0:1)
if(typ==0) hsyg=syg; Sigm=Sig; else hsyg=wsred; Sigm=Sigs; end
clear ws nbar x dbar;
for(hist_rownom=0:1)
if(hist_rownom==1)
lbsr=lbar/2; smax=max(hsyg); smin=min(hsyg);
ds=(smax-smin)/(lbar+1);
s=smin+ds/2; smax=smax-ds/2;
for(i=1:lbar) ws(i)=s; s=s+ds; end
[nbar,x]=hist(hsyg,lbar); %<=== nbar liczba wart. x w przedz.ds
nbar=nbar/ldan/ds; %<=== przelicz.czestosci na rozklad
else
lbmin=30; km=lbmin-2;
for(k=1:km) lbar=ldan/lbmin; if(lbar>=20) break; end;
lbmin=lbmin-1;
end
[nbar,x]=hist(hsyg,lbar); %<=== nbar l.wart. x w przedz.ds
k=0; s=0; xp=min(x);
for(i=1:lbar)
s=s+nbar(i);
if(s>=lbmin)
k=k+1; dbar(k)=s/ldan/(x(i)-xp); s=0;
ws(k)=(xp+x(i))/2; xp=x(i);
end
end
lbar=k; nbar=dbar; clear dbar;
end
for(i=1:lbar) fg(i)=exp(-((ws(i)/Sigm)^2)/2); end
fg=fg/Sigm/sqrt(2*pi);
subplot(2,2,hist_rownom+1+2*typ);
bar(ws,nbar);
hold on
plot(ws,fg,'r')
axis([min(ws)*1.05 max(ws)*1.05 0 max([max(fg) max(nbar)])]);
xlabel(sprintf('Rysunek %c',podpis(kk))); kk=kk+1;
clear x xs ws nbar fg;
end
end
end
35
Rys. 3 4 pokazują histogramy przykładowych danych losowych różnego typu.
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0
-1 0 1 -1 0 1
Rysunek a Rysunek b
2
2
1
1
0 0
-0.4 -0.2 0 0.2 -0.4 -0.2 0 0.2
Rysunek c Rysunek d
Rys.3. Histogramy 512 liczb losowych o rozkładzie równomiernym
a) histogram z nierównymi przedziałami; b) histogram z przedziałami równymi
c) histogram średnich z 20 próbek jak (a); d) histogram z przedziałami równymi jak (c)
36
0.2
0.2
0.15
0.15
0.1
0.1
0.05
0.05
0 0
-4 -2 0 2 4 -4 -2 0 2 4
Rysunek a Rysunek b
0.1
0.1
0.05
0.05
0 0
-5 0 5 -5 0 5
Rysunek c Rysunek d
Rys.4. Histogramy 512 liczb niejednorodnych losowych o rozkładzie równomiernym
a) histogram danych z dwiema różnymi wariancjami z nierównymi przedziałami;
b) jak (a) z przedziałami równymi; c) histogram danych z linioym wzrostem średniej
histogram nierównomierny jak (a); d) histogram danych jak (c)z przedziałami (b)
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0
-2 -1 0 1 -2 0 2
Rysunek a Rysunek b
2
2
1.5
1.5
1
1
0.5
0.5
0 0
-0.4 -0.2 0 0.2 0.4 -0.5 0 0.5
Rysunek c Rysunek d
Rys.5. Histogramy danych o rozkładzie normalnym jak na Rys.3.
37
0.2 0.15
0.1
0.1
0.05
0 0
-6 -4 -2 0 2 -5 0 5
Rysunek a Rysunek b
0.1 0.1
0.05 0.05
0 0
-6 -4 -2 0 2 4 -5 0 5
Rysunek c Rysunek d
Rys.6. Histogramy danych o rozkładzie normalnym jak na Rys.4.
38
Wyszukiwarka
Podobne podstrony:
Prezentacja ekonomia instytucjonalna na Moodle
model ekonometryczny zatrudnienie (13 stron)
Analiza ekonomiczna spółki Centrum Klima S A
Finanse Finanse zakładów ubezpieczeń Analiza sytuacji ekonom finansowa (50 str )
Wykład ekonomiczne podstawy
1 Wskaźniki techniczno ekonomiczne wiercenia otworuid049
Mysl Ekonomiczna i Polityczna 2 O Pietrewicz
Historia myli ekonomicznej wyklady
Ekonomia sektora publicznego 2010
Ekonomia Ebook Placet Ceny Tansferowe 1
więcej podobnych podstron