[vastaus aiempaan viestiin]
| Kirjoittaja: | Seppo Mustonen |
|---|---|
| Sähköposti: | - |
| Päiväys: | 5.4.2004 10:03 |
Tutkittuani erilaisia vaihtoehtoja liittää P-arvo REGDIAG:in
laskemaan Durbin-Watsonin testisuureeseen olen päätynyt seuraavaan
ratkaisuun, joka tulee mukaan versiosta 2.11 lähtien.
Katsotaan ensin, miten toimitaan käytännössä. Luon aluksi
näyteaineiston, jossa regressiomallin peräkkäiset residuaalit ovat
toisistaan riippuvia:
Aineiston DWKOE generointi:
Tehdään 4 muuttujan X1-X4 (teoreettinen) korrelaatiomatriisi R:
m=4 r=0.7
MAT R=CON(m,m,1-r)+IDN(m,m,r)
MAT RLABELS X TO R
MAT CLABELS X TO R
Generoidaan 100 havaintoa multinormaalijakaumasta N(0,R):
MNSIMUL R,*,DWKOE,100,0 / NEWSPACE=16,4
Lasketaan normaalisti jakautunut "selitysvirhemuuttuja" EPS:
VAR EPS=0.2*N.G(0,1,rand(2004)) TO DWKOE
Lasketaan selitettävä muuttuja Y:
VAR Y,RES TO DWKOE
Y1=1+X1+X2+X3+X4+EPS
Y=if(ORDER=1)then(Y1)else(Y1+0.3*EPS[-1])
RES=0
Teoreettiset regressiokertoimet ovat kaikki ykkösiä ja malliin liittyy
jäännöstermi EPS+0.3*EPS[-1],
jossa on jonkin verran "autokorrelaatiota".
.......................................................................
Regressiomallin laskenta:
MASK=XXXX-Y / Mallin valinta
REGDIAG DWKOE,CUR+1 / DWN=100000
Regression diagnostics on data DWKOE: N=100
Regressand Y # of regressors=5 (Constant term included)
Condition number of scaled X: k=1.67037
Variable Regr.coeff. Std.dev. t
Constant 0.9755829 0.0222096 43.926
X1 1.0097511 0.0254531 39.671
X2 0.9947873 0.0211140 47.115
X3 1.0059700 0.0248616 40.463
X4 1.0111429 0.0248550 40.682
Variance of regressand Y=5.847487900 df=99
Residual variance=0.045503393 df=95
R=0.9963 R^2=0.9925 Durbin-Watson=1.627 (P=0.0315) Autocorr>0
Yllä on aktivoitu REGDIAG täsmennyksellä DWN=100000.
Tällöin tavanomainen tulostus laajenee viimeisellä rivillä siten, että
saadaan myös Durbin-Watson-testisuureen P-arvo ja ilmoitus siitä, että
jäännöstermeissä on positiivista (1. kertaluvun) autokorrelaatiota.
P-arvon määrääminen tapahtuu Monte Carlo -menetelmällä (simuloimalla)
niin että koetoistoja on tässä tapauksessa 100000 (DWN=100000).
Omalla 1.6GHZ:in koneellani laskenta kesti noin 3.5 sekuntia.
Saatu P-arvo ei ole "tarkka", vaan se on luonteeltaan satunnainen,
keskivirheen ollessa s(P,DWN):=sqrt(P*(1-P)/DWN)
eli tässä tapauksessa s(0.0315,100000)=0.00055...
Todellinen P-arvo on siis joka tapauksessa suuruusluokkaa 0.03,
mikä tieto riittää hyvin tilastollisiin päätelmiin.
Teoreettista taustaa:
---------------------
Tarkastellaan lineaarista regressiomallia
Y = X*beta + eps,
missä Y (selitettävä) on n-vektori, X (selittäjät) on n*m-matriisi,
beta (regressiokertoimet) on m-vektori ja eps (selitysvirheet)
n-vektori, joka noudattaa jakaumaa N(0,sigma^2*I).
Kun merkitään
H=X*inv(X'*X)*X' ja M=I-H,
on
r=M*eps (mallin residuaalit)
ja Durbin-Watson testisuure
summa (r_i - r_i-1)^2 r'Wr eps' MWM eps
D = --------------------- = ---- = ------------ (koska MM=M)
summa (r_i)^2 r'r eps' M eps
missä W on yksinkertainen tridiagonaalinen n*n-matriisi.
Testisuureen D jakaumaa ei ole helppo käsitellä analyyttisesti,
sillä se riippuu olennaisesti matriisista X. Tästä johtuen
Durbin ja Watson aikoinaan (1950) päätyivät esittämään
testin kriittisille arvoille vain ala- ja ylärajoja, joita
nykyäänkin löytyy taulukoituna eri lähteissä.
Tuoreemmassa kirjallisuudessa esitetään kyllä keinoja D:n jakauman
kertymäfunktion laskemiseksi mutta silloin joudutaan esim. hankaliin,
suurten matriisien (kertalukua n*n) ominaisarvotehtäviin ja numeeriseen
integrointiin.
Nykyisillä, entistä nopeammilla koneilla on mielestäni käytännössä
järkevintä tyytyä määräämään testin P-arvo raa'alla simuloinnilla.
Menetellään näin:
Arvotaan eps, joka on N(0,I) ja lasketaan D em. neliömuotosuhteena.
Tätä toistetaan DWN kertaa ja lasketaan niiden toistojen suhteellinen
osuus, joissa D-arvo alittaa alkuperäisestä aineistosta lasketun arvon.
Tämä suhteellinen osuus lähestyy toistojen lisääntyessä testin P-arvoa.
Koska D on kahden neliömuodon (eps' MWM eps) ja (eps' M eps) suhde,
jossa kerroinmatriisit MWM ja M ovat suuria (siis n*n missä siis n
on havaintojen määrä), on katsottava tarkkaan, miten D esitetään
simuloitaessa. Suoraan laskettaessa kertolaskutoimituksia tulisi
(kun matriisit MWM ja M ovat valmiina) noin n^2 kpl.
Palataan sijoituksella r=M*eps (alkuperäiseen) muotoon
summa (r_i - r_i-1)^2 r'Wr
D = --------------------- = ----
summa (r_i)^2 r'r
ja käytetään hyväksi matriisin X singulaariarvohajotelmaa
X=U*D*V'
missä U on n*m-matriisi (U'U=I), D on m*m-lävistäjämatriisi ja
V m*m*-matriisi (V'V=I).
Tällöin "residuaaleille" r saadaan lauseke
r = M*eps = (I-H)*eps = (I-X*inv(X'X)*X')*eps
= (I-U*U')*eps
= eps - U*U'*eps.
Siten D:n laskentaan tarvittava kertolaskujen lukumäärä on enää
suuruusluokkaa 2*m*n, mikä on suurilla havaintomäärillä n olennaisesti
pienempi kuin n^2.
Tämä on Survon yhteydessä muutenkin luonnollinen tie, sillä REGDIAG
käyttää em. hajotelmaa X=U*D*V' jo itse regresiolaskelmissa. Se on
tavanomaiseen normaaliyhtälöiden ratkaisumenetelmään verrattuna
tarkempi ja turvallisempi menettely.
* * *
Ne, joita kiinnostaa teoreettisemmat lähestymistavat, voivat
tutkiskella artikkeleita (sain nämä vihjeet Pentti Saikkoselta):
Shively, T.S., Ansley, C.F. & Kohn, R. (1990). Fast evaluation of the
distribution of the Durbin-Watson and other invariant test statistics in
time series regression. Journal of the American Statistical Association
85, 676-685.
Ansley, C.F., Kohn, R. & Shively, T.S. (1992). Computing p-values for
the generalized Durbin-Watson and other invariant test statistics.
Journal of Econometrics 54, 277-300.
Ideana P-arvojen laskemisessa näyttää olevan testisuureen jakauman
karakteristisen funktion kääntäminen numeerisesti.
* * *
Olen kuitenkin jo pitkän aikaa (jo 1980-luvun alusta lähtien) suosinut
Survossa satunnaistettuja testejä. Sellaisesta on edellä esittämässäni
laskentatavassakin kysymys. Näin tullaan usein toimeen huomattavasti
vähemmillä olettamuksilla. Vältetään myös liikaa teoreettista
painolastia ja laskuteknisiä ongelmia.
Uskon, että tulevaisuus koituu yhä vahvemmin satunnaistettujen
menettelyjen eduksi.
-Seppo
| Vastaukset: |
|---|
Survo-keskustelupalstan (2001-2013) viestit arkistoitiin aika ajoin sukrolla, joka automaattisesti rakensi viesteistä (yli 1600 kpl) HTML-muotoisen sivukokonaisuuden. Vuoden 2013 alusta Survo-keskustelua on jatkettu entistäkin aktiivisemmin osoitteessa forum.survo.fi. Tervetuloa mukaan!