Дискрет сонголтын загвар - Discrete Choice Models
Д. Хашбаатар
Бусад нэршлүүд: Limited dependent variable, qualitative variable, categorical variable
- 1 Binary Choice (Logit, Probit)
- 2 Multinomial Choice (Multinomial Logit, Probit)
- 3 Count Data (Poisson Regression)
- 4 Truncation, Censor, Duration et al (Tobit, Heckit, Cox, Weibull)
1 Binary Choice Model
1.1 хоёр сонголтын жишээ
- сommuting: нийтийн тээвэр / хувийн машин
- сонгууль: АН / МАН
- үүрэн телефон: Мобиком / Скайтел
- зээл: зөвшөөрсөн / татгалзсан
Ийм төрлийн шинжилгээ ихэвчлэн хувь хүн, өрх, компанийн түвшний дата дээр суурилсан байдаг учир микроэконометрик гэж ангилан үзэх нь бий.
1.2 Linear Probability Model
Өрхүүд машин эзэмших нь гэр бүлийн орлогоос хэрхэн хамаарч буйг тайлбарлах загвар зохиоё. Тайлбарлагч хувьсагч нь гэр бүлийн орлого \(x_i\). Машин худалдан авсан эсэхийг илтгэгч хувьсагч
\[ y_i= \begin{cases} 1& i\text{-р өрх машинтай}\\ 0& i\text{-р өрх машингүй} \end{cases} \]
Үүнийг шугаман регрессээр тайлбарлах гэж оролдвол \[y_i=\beta_0+\beta_1 x_i+u_i\] загварыг үнэлэн коэффициентуудыг олж болох боловч илэрхийлэх утга нь энгийн шугаман регрессээс өөр байна. Тухайлбал,
- \(\beta_1=\Delta y/\Delta x\)?
- \(\beta_0+\beta_1 x\) шулуун ямар утгатай вэ?
- \(\hat{y}\)-н утга? \(\hat{y}=0.3\)-г юу гэж тайлбарлах вэ?
Шугаман регрессийн хувьд \[E(y|x)=\beta_0+\beta_1 x\] байдаг.
\(y\) нь 2 утгатай тул \[E(y|x)=P(y=1|x)1+P(y=0|x)0=P(y=1|x)\] учир
\[P(y=1|x)=\beta_0+\beta_1 x_i\]
Өөрөөр хэлбэл \(y\) нь бинар хувьсагч бол
- \(\beta_0+\beta_1 x\) шулуун нь \(y=1\) байх магадлалыг илэрхийлнэ.
- \(\hat{y}\) нь таамаглагдсан магадлал (predicted probability)
Алдаа нь heteroskedastic учир heteroskedastic-robust вариацын үнэлгээ ашиглах нь зүйтэй.
Шугаман регрессийн сул талууд:
- \(\hat{y}>1\) юмуу \(\hat{y}<0\)
1.3 Ложит ба Пробит (Logit / Probit)
\[E(y)=G(\beta_0 +\beta_1 x)\] энд \(0\leq G\leq 1\) байх зааглагдсан функц гэж сонгож авбал \(\hat{y}\) үргэлж магадлал илэрхийлнэ.
- Probit: \[G(x)=\Phi(x)\]
- Logit: \[G(x)=\frac{\exp(x)}{1+\exp(x)}\]
Ложит болон пробит функцууд хоорондоо төстэй учир үр дүн нь ч ойролцоо байдаг.
\[\frac{\partial y}{\partial x}=\frac{\partial \Phi(\beta_0+\beta_1x)}{\partial x}=\phi(\beta_0+\beta_1x)\beta_1\]
1.3.1 Latent variable interpretation
Далд хувьсагчийн тайлбар: Ханамж максимумчлах бодлогоос ложит/пробит моделийг гарган авч болдог.
Хэрэглэгчийн сонголтын олонлог 2 элементтэй, \(y=\{1, 0\}\). Хэрэглэгчийн ханамжийг эконометрикч ажиглаж чадахгүй боловч WARP (weak axiom of revealed preferences) ёсоор \(y=1\) бол \(u(1)-u(0)>0\) байх ёстой. Ханамжийн ялгааг эконометрикч \[y^*=x'\beta+\varepsilon\] гэж загварчилбал \[P(y=1)=P(y^*>0)=P(-\varepsilon<x'\beta).\] \(\varepsilon\) нь симметрик тархалттай бол дээрхи хэмжигдэхүүн \(F(x'\beta)\) болно. Нормаль тархалттай гэж үзвэл пробит модел гарна.
1.4 Үнэлгээ
Likelihood function \[L(\beta)=\prod_{i=1}^nP(y=1|x_i)^{y_i}P(y_i=0|x_i)^{1-y_i}\]
\(P(y=1|x_i)=F(x_i'\beta)\) гэж орлуулан логарифм авбал (тэмдэглэгээ багахан өөрчлөгдснийг анхаар!)
\[\log L(\beta)=\sum_{i=1}^ny_i\log P(y=1|x_i)+\sum_{i=1}^n(1-y_i)\log P(y_i=0|x_i)\]
FOC for maximization
\[\frac{\partial \log L(\beta)}{\partial \beta}= \sum_{i=1}^n\left[\frac{y_i-F(x'\beta)}{F(x'\beta)(1-F(x'\beta))}f(x'\beta)\right]f(x'\beta)=0\]
OLS-тэй адил аналитик шийд байхгүй учир дөхөх алгоритмаар шийдийг олно.
2-р эрэмбийн уламжлалууд сөрөг тодорхойлогдсон матриц үүсгэх учир, зорилтын функц маань (глобаль) гүдгэр болохыг шалгаж болно. Иймд цор ганц максимумтай бөгөөд Ньютон-Рапсон зэрэг дөхөх алгоритм зайлшгүй нийлэх болно.
Хамгийн их үнэний утгын үнэлгээ (MLE) нь асимптот нормаль учир OLSтэй адилаар t-тест, F-тест ашиглан таамаглал шалгаж болно.
1.5 Таарамж (Goodness of Fit)
Регрессанд \(y\) нь дискрет (or чанарын) хувьсагч тул \(R^2\)-г шугаман регрессийн адил тодорхойлох боломжгүй.
pseudo R^2
\(\log L_1\): maximum log likelihood, \(\log L_0\): log likelihood when all parameters are 0, i.e Bernoulli maximum likelihood.
- \[pseudoR^2 = 1 − \frac{1}{1 + 2(\log L_1 − \log L_0)/n}\]
- \[McFaddenR^2 = 1 − \frac{\log L_1}{\log L_0}\]
fraction correctly predicted uses following rule: \(y=1\) if \(F(x'\beta)>1/2\) gevel \(\hat{y}=1\) if \(x'\hat\beta>0\)
\[R_p^2=1-\frac{wr_1}{wr_0}\]
1.6 R Implementation
Статистикт \[E(y|x)=g(x'\beta)\] моделийг GLM(Generalized Linear Model) гэж нэрлэдэг. Эконометрикт ихэвчлэн шугаман бус регресс (nonlinear regression) нэрээр тааралдана. \(g\)-нь \(y\)-ийн(regressand) дунджыг \(x\)-ийн(regressor) шугаман функцтэй холбож буй учир линк функц (link function) гэж нэрлэдэг. Binary choice модель маань GLM-ийн зөвхөн нэг хэсэг.
R программ дээр үнэлгээг glm коммандаар хийнэ. Үр дүн нь шугаман регрэссийн lm object -тэй төстэй list байдлаар гарна. family (binomial, gaussian, poisson,inverse.gaussian гэх мэт), мөн түүн дотроо link( probit, logit,log,identity гэх мэт)-ийг сонгох сонголттой. Жишээ нь пробит моделийг доорх байдлаар үнэлнэ.
probit_est <- glm(y ~ x, data = df, family = binomial(link = "probit"))
summary(probit_est)1.7 Жишээ
АНУ-д (зээл олгоход) арьс өнгөөр ялгаварлахыг хуулиар хориглосон байдаг. Гэтэл 1990 оны Бостонны Массачусетсийн нийт зээл хүсэгчдийн өгөгдөлөөс харахад хар арьстаны 28% нь зээл авч чадаагүй бол цагаан арьстаны 9% нь л зээл авч чадаагүй байна.
Бостон HMDA (Home Mortgage Disclosure Act ):2Munnell, A. H., G. M. B. Tootell, L. E. Browne, and B. J. McEneaney (1996): Mortgage Lending in Boston: Interpreting HMDA Data, American Economic Review, 86, 25–53. Бостоны ипотекийн зээлийн дата.
library(readxl)
suppressMessages(library(dplyr))
library(magrittr)
library(ggplot2)
hmda <- read_excel("datasets/hmda_sw.xlsx")
hmda %<>%
mutate(black = ifelse(s13==3, 1, 0),
deny = ifelse(s7==3, 1, 0),
PI = s46/100)
hmda %>%
group_by(black, deny) %>%
summarise(too=n()) ## Source: local data frame [4 x 3]
## Groups: black [?]
##
## black deny too
## (dbl) (dbl) (int)
## 1 0 0 1852
## 2 0 1 189
## 3 1 0 243
## 4 1 1 96f1 <- lm(deny ~ PI, data = hmda)
summary(f1)##
## Call:
## lm(formula = deny ~ PI, data = hmda)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73070 -0.13736 -0.11322 -0.07097 1.05577
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.07991 0.02116 -3.777 0.000163 ***
## PI 0.60353 0.06084 9.920 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3183 on 2378 degrees of freedom
## Multiple R-squared: 0.03974, Adjusted R-squared: 0.03933
## F-statistic: 98.41 on 1 and 2378 DF, p-value: < 2.2e-16
\[\hat{deny}=-0.08+0.604PI\]
f2 <- lm(deny ~ PI+black, data = hmda)
summary(f2)##
## Call:
## lm(formula = deny ~ PI + black, data = hmda)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62526 -0.11772 -0.09293 -0.05488 1.06815
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.09051 0.02079 -4.354 1.39e-05 ***
## PI 0.55919 0.05987 9.340 < 2e-16 ***
## black 0.17743 0.01837 9.659 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3123 on 2377 degrees of freedom
## Multiple R-squared: 0.076, Adjusted R-squared: 0.07523
## F-statistic: 97.76 on 2 and 2377 DF, p-value: < 2.2e-16\[\hat{deny}=-0.09+0.559PI+0.177black\]
f3 <- glm(deny ~ PI, data=hmda, family=binomial(link="probit"))
f3logit <- glm(deny ~ PI, data=hmda, family=binomial(link="logit"))
attributes(f3)## $names
## [1] "coefficients" "residuals" "fitted.values"
## [4] "effects" "R" "rank"
## [7] "qr" "family" "linear.predictors"
## [10] "deviance" "aic" "null.deviance"
## [13] "iter" "weights" "prior.weights"
## [16] "df.residual" "df.null" "y"
## [19] "converged" "boundary" "model"
## [22] "call" "formula" "terms"
## [25] "data" "offset" "control"
## [28] "method" "contrasts" "xlevels"
##
## $class
## [1] "glm" "lm"print(object.size(f3), units = "MB")## 2.4 Mbbroom::tidy(f3)## term estimate std.error statistic p.value
## 1 (Intercept) -2.194145 0.1377614 -15.927140 4.107187e-57
## 2 PI 2.967867 0.3857570 7.693617 1.430332e-14broom::tidy(f3logit)## term estimate std.error statistic p.value
## 1 (Intercept) -4.028432 0.2685756 -14.999250 7.425337e-51
## 2 PI 5.884498 0.7335988 8.021411 1.045371e-15
\[\widehat{P(deny=1|PI)}=\Phi(-2.19+2.97PI)\] \[\widehat{P(deny=1|PI)}=L(-4.02+5.88PI)\]
newdata = data.frame(PI=c(0.3, 0.4))
predict(f3, newdata, type="response")## 1 2
## 0.09615344 0.15696777predict(f3logit,newdata, type="response")## 1 2
## 0.09422692 0.15780744
f4 <- glm(deny ~ PI+black, data=hmda, family=binomial(link="probit"))
broom::tidy(f4)## term estimate std.error statistic p.value
## 1 (Intercept) -2.2587870 0.13669056 -16.524821 2.431605e-61
## 2 PI 2.7417794 0.38046928 7.206310 5.748862e-13
## 3 black 0.7081554 0.08335249 8.495911 1.963866e-171.7.1 Descriptive Statistics of Boston HMDA
hmda %<>%
mutate(hse_inc = s45/100,
loan_val = s6/s50,
ccred = s43,
mcred = s42,
pubrec = 1*(s44>0),
denpmi = 1*(s53==1),
selfemp = 1*(s27a==1),
married = 1*(s23a=="M"),
single = 1*(married==0),
hischl = 1*(school>=12),
probunmp = uria,
condo = 1*(s51 == 1))
means <- hmda %>%
summarise_each(funs(mean), PI, hse_inc, loan_val, ccred, mcred, pubrec, denpmi, selfemp,
single, hischl, probunmp, condo, black, deny)
reshape2::melt(means)## No id variables; using all as measure variables## variable value
## 1 PI 0.33081357
## 2 hse_inc 0.25534612
## 3 loan_val 0.73777591
## 4 ccred 2.11638655
## 5 mcred 1.72100840
## 6 pubrec 0.07352941
## 7 denpmi 0.02016807
## 8 selfemp 0.11638655
## 9 single 0.39327731
## 10 hischl 0.98361345
## 11 probunmp 3.77449585
## 12 condo 0.28823529
## 13 black 0.14243697
## 14 deny 0.119747902 Нэмэлт материал
Kenneth Train “Discrete Choice Methods with Simulation”
Imbens/Wooldridge “Lecture 11: Discrete Choice Models” NBER Summer Institute Econometric Lectures 2007