Bayesiánská analýza - cvičení 3
Toto cvičení je založeno na znalosti čtvrté kapitoly z učebnice Koop (2003): Bayesian econometrice, případně na odpovídající kapitole podkladového učebního textu Bayesiánská analýza.
Co bude náplní cvičení?
^ Odhad a posteriorní analýza normálního lineárního regresního modelu s nezávislou normální-gama apriorní hustotou a s omezeními ve tvaru nerovnosti.
^ Osvojení si Gibbsova vzorkovače. Zadání příkladů
1. Gibbsův vzorkovač a jeho vlastnosti: Předpokládejme elementární příklad modelu, kdy při odhadu jeho parametrů získáváme posteriorní hustotu odpovídající dvourozměrnému normálnímu rozdělení (to, že máme známou posteriorní hustotu bude sloužit k tomu, že výsledky simulátorů budeme schopni porovnat s výsledky analytickými):
kde \p\ < 1 je známá posteriorní korelace mezi parametry 6\ a 62-
(a) Začněte si vytvářet skript, kdy si nejprve zadefinujete počet generovaných vzorků S (např. S = 10000), vytvořte vektor středních hodnot p jakožto nulový vektor (později ho můžete změnit), korelační koeficient p a kovarianční matici X, kdy na diagonále budou jedničky a mimo ní koeficienty p (díky tomu odpovídá kovarianční matice přímo i korelační matici).
(b) Vytvořte si vlastní funkci, která využívá Monte Carlo integraci k výpočtu posteriorní střední hodnoty, směrodatné odchylky parametrů 6\ a 62 a numerické standardní chyby (NSE) pro střední hodnotu a rozptyl odhadu parametrů, popřípadě využijte dodanou funkci MC_int. m. V rámci této funkce můžete využít generátor náhodných čísel z vícerozměrného rozdělení buď přímo dostupný v rámci statistického tool-boxu Matlabu (funkce mvnrnd.m) nebo funkci LeSageho ekonometrického toolboxu norm_rnd.m. Poznámka: Je třeba si uvědomit, že Monte Carlo integrace nevyužívá nic jiného než zákon velkých čísel, který říká, že výběrový průměr nějaké funkce parametrů bude konvergovat ke střední hodnotě této funkce parametrů pro rostoucí velikost vygenerovaných vzorků rozdělení, z něhož tyto parametry pocházejí. NSE
je definováno jako y    , kde a'1 je rozptyl funkce parametrů (nahrazován odhadem), který nás zajímá.
(c) Vytvořte funkci, která budeý využívat Gibbsův vzorkovač k výpočtu střední hodnoty a směrodatné odchylky parametrů 6\ a 62 (využijte vlastností vícerozměrného normálního rozdělení pro výpočet odpovídajících podmíněných hustot pravděpodobnosti).
i. Vyjděte z funkce sdružené hustoty normálního rozdělení (dvojrozměrné), se střední hodnotou p = (/líi, P2)' a kovarianční maticí X definovanou v úvodu příkladu:
ii. Odvoďte podmíněné hustoty pro 6\ #2 a #2 |#i
iii. Pokud jste šikovnější, odvoďte podmíněné hustoty pro obecné fc-rozměrné normální rozdělení s arbit-rárním dělením vektoru středních hodnot a kovarianční matice (zde je dobré využít teorém o inverzi a determinantu dělené matice z přílohy učebního textu či Koopovy učebnice). Důkaz inverze dělené matice lze nalézt např. zde a odvození podmíněných hustot zde.
1
£|-*exp --(e-py^-He-p)
{2tt)Í
1
iv. Výsledné funkce podmíněných hustot pravděpodobnosti jsou 6\#2 ~ -^(M(i|2)j £(i|2)) a ®2\Qi ~ iV(M(2|i),S(2|i)):
M(i|2) = Mi +P* (#2 -M2)
S(l|2) = 1-P2
M(2|i) = M2 + P* (0i - Mi)
£(2|i) = 1-P2
(d) Nastavte p = 0 pro porovnání výsledků z části (a) a (b). Kolik replikací je nutných k odhadu středních hodnot a směrodatných odchylek parametrů 6\ a 62 s přesností na dvě desetinná místa?
(e) Zopakujte část (c) pro p = 0.5, 0.9, 0.95, 0.99, 0.999. Jak velikost korelace ovlivní výkonnost Gibbsova vzorkovače? Pro srovnání vykreslete průběh Gibbsova vzorkovače pro prvních 50 a 1000 iterací a rovněž i výslednou sekvenci vzorků po odstranění počátečních Sq vzorků, a to pro odlehlé počáteční hodnoty parametrů.
2. Soubor data_cocaine . mat obsahuje 56 pozorování proměnných vztahujících se k prodeji kokainu v severovýchodní Kalifornii v období 1984-1991. Data jsou podmnožinou dat použitých ve studii Culkins, J.P. a Pad-man, R. (1993): „Quantity Discounts and Quality Prémia for Illicit Drugs," Journal of the American Statistical Association, 88, 748-757. Proměnné jsou načteny v rámci výchozího skriptu cocaine . m a mají následující význam:
• price = cena za gram kokainu v rámci dané transakce;
• quant = počet gramů kokainu prodaných v dané transakci;
• qual = kvalita kokainu vyjádřená jako procento čistoty;
• trend = časová proměnná s hodnotami od 1984=1 až po 1991=8.
Předpokládejme regresní model
price = /3q + /3iquant + faqual + fi^trend + e.
(a) Jaká znaménka koeficientů byste očekávali u parametrů Pi, P2 a fiz?
(b) Odhadněte daný model (předpokládáme, že se jedná o NLRM s nezávislou normální gama apriorní hustotou). Zvolte si vhodné hyperparametry dle vašich zkušeností. Jsou znaménka parametrů v souladu s vašim očekáváním?
(c) Říká se, že čím větší objem obchodů, tím větší riziko, že vás dostihne ruka zákona. Prodejci tak jsou ochotni akceptovat nižší cenu, pokud prodávají větší množství. Pokuste se testovat tuto hypotézu.
(d) Ověřte hypotézu, že kvalita kokainu nemá vliv na jeho cenu.
(e) Jaká je průměrná roční změna ceny kokainu? Zamyslete se nad tím, proč by se měla cena takto měnit.
(f) Odhadněte výchozí model se zahrnutím o předpokladech na znaménka parametrů dle vašich očekávání a otestujte znovu vybrané hypotézy.
2
3. Každé ráno mezi 6:30 a 8:00 opouští Bill Melbournské předměstí Carnegie, aby se dostal do práce na University of Melbourne. Čas, který Bili stráví cestou do práce, Ume, závisí na času odjezdu, depart, počtu červených světel na semaforech, reds a počtu vlaků, kvůli kterým musí čekat na Murrumbeenském přejezdu, trains. Pozorování těchto proměných je celkem získáno za 231 pracovních dní v roce 2006 a jsou obsahem souboru data_commute . mat, resp. jsou jíž načtena v rámci skriptu commute . m. Proměnná Ume je měřena v minutách, depart je počet minut po 6:30, které uplynou než Bili vyrazí z domu.
(a) Odhadněte rovnici (v kontextu NLRM s nezávislou normální gama apriorní hustotou)
Ume = ßo + ßidepart + ß'ireds + ß^trains + e.
(b) Jaká znaménka koeficientů byste očekávali u parametrů ß\, ßi a ßs?
(c) Otestujte hypotézu, že každé červené světlo zpozdí Billa nejméně o 2 minuty.
(d) Testujte hypotézu, že čas odjezdu nemá vliv na čas strávený cestováním.
(e) Otestujte hypotézu, čas cestování navíc díky čekání na jednom semaforu je stejný jako čas čekání průjezdu jednoho vlaku.
(f) Odhadněte výchozí model se zahrnutím o předpokladech na znaménka parametrů dle vašich očekávání a otestujte znovu vybrané hypotézy.
3