/* Na podstawie przykładów Cameron,Trivedi (2009) Microeconometric Using Stata" */
* Załadowanie danych - są zamieszczone na stronie Stata Corp.
net from http://www.stata-press.com/data/mus
use mus14data.dta
* definiowanie interakcji
drop age2 agefem agechr agewhi
global xlist age hstatusg hhincome educyear married hisp
generare linc=ln(hhinc)
global extralist linc female white chronic adl sretire

summarize ins retire $xlist $extralist

*Model logitowy
logit ins retire $xlist
estimates store blogit

*Model probitowy
qui probit ins retire $xlist
estimates store bprobit

*Liniowy model prawdopodobieństwa
qui reg ins retire $xlist
estimates store bols

*Logit z macierzą odporną
qui logit ins retire $xlist, vce(robust)
estimates store blogitr

*Probit z macierzą odporną
qui probit ins retire $xlist, vce(robust)
estimates store bprobitr

*Regresja z macierzą odporną
qui reg ins retire $xlist, vce(robust)
estimates store bolsr

*Podsumowanie i porównaie modeli
estimates table blogit blogitr bprobit bprobitr bols bolsr, t stats(N llo) b(%7.3f) stfmt(%8.2f) 

* Test Walda (istotność interakcji)
generate age2 = age*age
generate agefem = age*female
generate agechr = age*chronic
generate agewhi = age*white

global intlist age2 agefem agechr agewhi

quietly logit ins retire $xlist $intlist
test $intlist

*Test LR (istotność interakcji)
quietly logit ins retire $xlist $intlist
estimates store B

quietly logit ins retire $xlist 
lrtest B

/* Test poprawności formy funkcyjnej */
quietly logit ins retire $xlist
linktest

hetprob ins retire $xlist, het(chronic) //Heteroscedastic Probit

/* Miary dopasowania */

quietly logit ins retire $xlist
estat gof, group(4)
estat gof, group(10)

estat classification

/* Predykcja */
* calculate predicted probabilities
quietly logit ins hhincome
predict plogit, pr

quietly probit ins hhincome
predict pprobit, pr

quietly reg ins hhincome
predict pols, xb

sort hhincome
graph twoway (scatter ins hhincome,msize(vsmall) jitter(3)) /*
*/ (line plogit hhincome, clstyle(p1) ) /*
*/ (line pprobit hhincome, clstyle(p2) ) /*
*/ (line pols hhincome, clstyle(p3) ) /*
*/ scale(1.2) plotregion(style(none)) /*
*/ title("Predicted probabilities across models") /*
*/ xtitle("HHINCOME (hhincome)", size(medlarge)) xscale(titlegap(*5)) /*
*/ ytitle("Predicted probability", size(medlarge)) yscale(titlegap(*5)) /*
*/ legend(pos(1)) ring(0) col(1)) legend(size(small))  /*
*/ legend(label(1 "Actual Data (jittered)") label(2 "Logit") label(3 "Probit") label(4 "OLS"))

/* efekty krańcowe */
quietly logit ins retire $xlist
mfx
mfx, at(1 75 1 35 12 1 1)
/* przecietna wielkość pochodnej */
net serach mergeff
margeff
/* zmiana efektów krańcowych wywołana jednostkową zmianą regresora */
prchange

/* model ze zmiennymi endogenicznymi */
global xlist2 female age age2 educyear married hisp white chronic adl hstatusg
probit ins linc $xlist2, vce(robust)

global ivlist retire sretire
ivprobit ins $xlist2 (linc=$ivlist),vce(robust)

/* Intepretacja interakcji */
sysuse nlsw88

generate byte high_occ = occupation < 3 if occupation < .
generate byte black = race == 2 if race < .
drop if race == 3
generate byte baseline = 1
logit high_occ black##collgrad baseline, or noconstant nolog


margins, over(black collgrad) expression(exp(xb())) post