setwd("C:/Users/marcogeovanni/Desktop/Modern Guide To Econometrics/Chapter 7")
library(readstata13)
CreditRatings <-read.dta13("credit.dta")
attach(CreditRatings)
summary(CreditRatings)
## booklev ebit invgrade logsales
## Min. :0.0000 Min. :-0.38417 Min. :0.0000 Min. : 1.100
## 1st Qu.:0.1698 1st Qu.: 0.05120 1st Qu.:0.0000 1st Qu.: 6.970
## Median :0.2640 Median : 0.09044 Median :0.0000 Median : 7.884
## Mean :0.2932 Mean : 0.09389 Mean :0.4723 Mean : 7.996
## 3rd Qu.:0.3877 3rd Qu.: 0.13599 3rd Qu.:1.0000 3rd Qu.: 8.950
## Max. :0.9992 Max. : 0.65151 Max. :1.0000 Max. :12.701
## marklev rating reta wka
## Min. :0.0000 Min. :1.000 Min. :-0.99589 Min. :-0.41208
## 1st Qu.:0.1091 1st Qu.:3.000 1st Qu.: 0.01108 1st Qu.: 0.02905
## Median :0.2112 Median :3.000 Median : 0.18048 Median : 0.12284
## Mean :0.2547 Mean :3.499 Mean : 0.15699 Mean : 0.14041
## 3rd Qu.:0.3477 3rd Qu.:4.000 3rd Qu.: 0.35084 3rd Qu.: 0.23462
## Max. :0.9649 Max. :7.000 Max. : 0.97992 Max. : 0.74802
library(MASS)
## Warning: package 'MASS' was built under R version 3.2.3
library(erer)
## Warning: package 'erer' was built under R version 3.2.3
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 3.2.3
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
Logit <- glm(invgrade~booklev + ebit + logsales + reta + wka, family = binomial(link= "logit"))
summary(Logit)
##
## Call:
## glm(formula = invgrade ~ booklev + ebit + logsales + reta + wka,
## family = binomial(link = "logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5840 -0.5411 -0.0491 0.5286 3.4210
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.21432 0.86685 -9.476 < 2e-16 ***
## booklev -4.42727 0.77142 -5.739 9.52e-09 ***
## ebit 4.35474 1.43992 3.024 0.00249 **
## logsales 1.08159 0.09568 11.304 < 2e-16 ***
## reta 4.11611 0.48851 8.426 < 2e-16 ***
## wka -4.01249 0.74791 -5.365 8.10e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1273.95 on 920 degrees of freedom
## Residual deviance: 682.16 on 915 degrees of freedom
## AIC: 694.16
##
## Number of Fisher Scoring iterations: 6
LogisticDensityFunction <- mean(dlogis(predict(Logit, type= "link")))
LogisticDensityFunction
## [1] 0.1184878
LogisticDensityFunction*coef(Logit)
## (Intercept) booklev ebit logsales reta wka
## -0.9732966 -0.5245769 0.5159829 0.1281555 0.4877085 -0.4754313
Logit0 <- update(Logit, formula = invgrade~1 )
Mcfadden <-1-(logLik(Logit)/logLik(Logit0))
print(Mcfadden)
## 'log Lik.' 0.464536 (df=6)
plot(rating)
ratingfactor <-factor(rating)
OrderedLogit <- polr(ratingfactor~booklev + ebit + logsales + reta + wka)
summary(OrderedLogit)
##
## Re-fitting to get Hessian
## Call:
## polr(formula = ratingfactor ~ booklev + ebit + logsales + reta +
## wka)
##
## Coefficients:
## Value Std. Error t value
## booklev -2.7517 0.47730 -5.765
## ebit 4.7313 0.94476 5.008
## logsales 0.9406 0.05884 15.984
## reta 3.5600 0.30229 11.777
## wka -2.5802 0.48263 -5.346
##
## Intercepts:
## Value Std. Error t value
## 1|2 -0.3696 0.6332 -0.5837
## 2|3 4.8805 0.5208 9.3709
## 3|4 7.6259 0.5512 13.8361
## 4|5 9.8849 0.5916 16.7078
## 5|6 12.8829 0.6734 19.1315
## 6|7 14.7825 0.7841 18.8527
##
## Residual Deviance: 1930.614
## AIC: 1952.614
OrderedLogit0 <- update(OrderedLogit, formula = ratingfactor~1)
print(OrderedLogit)
## Call:
## polr(formula = ratingfactor ~ booklev + ebit + logsales + reta +
## wka)
##
## Coefficients:
## booklev ebit logsales reta wka
## -2.7516633 4.7313154 0.9405506 3.5599921 -2.5801607
##
## Intercepts:
## 1|2 2|3 3|4 4|5 5|6 6|7
## -0.3696357 4.8805454 7.6258969 9.8849375 12.8828577 14.7825479
##
## Residual Deviance: 1930.614
## AIC: 1952.614
ocME(OrderedLogit)
##
## Re-fitting to get Hessian
## effect.1 effect.2 effect.3 effect.4 effect.5 effect.6 effect.7
## booklev 0.001 0.196 0.478 -0.487 -0.177 -0.009 -0.002
## ebit -0.002 -0.337 -0.822 0.838 0.305 0.016 0.003
## logsales 0.000 -0.067 -0.163 0.167 0.061 0.003 0.001
## reta -0.002 -0.254 -0.619 0.631 0.229 0.012 0.002
## wka 0.001 0.184 0.448 -0.457 -0.166 -0.009 -0.002
McfaddenOrdered <-1-logLik(OrderedLogit)/logLik(OrderedLogit0)
print(McfaddenOrdered)
## 'log Lik.' 0.3088874 (df=11)