#HANDS-ON INTERMEDIATE ECONOMETRICS USING R 
# Second Edition, 2022, World Scientific Publishers
#by Hrishikesh D Vinod 
#Following snippets are from Chapter 7
#If you use them, please give credit to Prof. Vinod



#R7.1.1 LPM for Boston mortgage denial and family indebtedness
library(Ecdat); data(Hdma)
dat=data.frame(Hdma)
summary(dat)
names(dat)
n=length(dat$deny)
X=dat$dir  #debt payments to total income ratio
Y=rep(0,n) #initialize space to store Y
for (i in 1:n) {
  if (dat$deny[i]=="yes")  Y[i]=1} #Y=1 if mortgage is denied
# Y=as.numeric(dat$deny)-1 #avoids the loop above
reg1=lm(Y~X) # intercept =-0.07996, slope= 0.60353
f=fitted(reg1)  #collect fitted values of the regression
cor(f,Y)^2# 0.03973527 is the R-square reported by OLS
#function lm (squared correlation bet Y and fitted Y)
length(f[f<0]) #count 39 fitted values are negative
length(f[f>1])# 1 fitted value exceeds 1, so cannot be a
#probability
library(car)
ncv.test(reg1) #tests Null of constant error variance against
#changing heteroscedastic variance
#Chisquare = 116.4283    Df = 1     p = 0
library(lmtest) # perform Breusch-Pagan test
bptest(reg1) #BP = 4.4939, df = 2, p-value = 0.1057
hmctest(reg1, order.by=X) #Harrison-Mccabe test p-value
#< 2.2e-16 means we reject homoscedasticity


#R7.1.2 Probit for Boston mortgage denial and family
#indebtedness
# Assume that #R7.1.1 is already run and in the memory of R
pro=glm(Y~X, family=binomial(link="probit"))
summary(pro)
fp=fitted(pro)  #collect fitted probit regr. values
length(fp[fp<0]) #No. of negative fitted values=0
length(fp[fp>1]) #No. of fitted values larger than 1=0
cor(fp,Y)^2 # 0.05152176
dydxp=pro$coef[2]*dnorm(fp)#marginal effect of the slope
md=mean(dydxp); md #marginal effect on prob eval at mean
# is 17% larger than unit
md*(sd(X)/10)#marinal effect of unit=sd/10 is 17%
#larger than unit



#R7.1.3 Logit for Boston mortgage denial and family
#indebtedness
# Assume that #R7.1.1 is already run and in the memory of R
logi=glm(Y~X, family=binomial(link="logit"))
summary(logi)
fL=fitted(logi)  #collect fitted values of the logit
#regression
length(fL[fL<0])#count No. of negative fitted values
length(fL[fL>1]) #count fitted values larger than 1
cor(fL,Y)^2 # 0.05152176 goodness of fit like R-square
dydxL=logi$coef[2]*fL*(1-fL)#marginal effect of the slope
mean(dydxL) #mean of marginal effects= 0.5907.
avOddsR=mean(exp(fL))  #average odds ratio in the data
chg=avOddsR-1 #change in odds ratio
dydx.OddsR=chg*logi$coef[2]
mean(dydx.OddsR) #marg. effect on odds ratio (subtracting 1)



#R7.1.4 Example:grade improvement by new instruction method
#source Spector and Mazzeo, J Econ Education vol.11, 1980,
#pp37-44 also available at link li below.
li="https://faculty.fordham.edu/vinod/logit.txt"
econgrade=read.table(file=li,header=TRUE)
head(econgrade,2); tail(econgrade,2)
#OBS, gpa (grade point average), ecoscore (score in econ test)
#particip (binary variable 1=participated in new method)
# gradeincr (binary variable 1=yes grade did increase)
lgt=glm(gradeincr~gpa+ecoscore+particip, family=binomial
        (link="logit"), data=econgrade )
summary(lgt)  #results match with Greeene text
#(Ed.5)page 677.
pbt=glm(gradeincr~gpa+ecoscore+particip, family=binomial
        (link="probit"), data=econgrade)
summary(pbt)


#R7.1.5  Typical set up of logit models in R
#This will not work as is, user must supply y, X data
#glm(y~x1+x2+x3, family=binomial("logit"))
#where y is often a matrix with 2 columns for success vs.
#failure, or dead vs. alive.
#Other choices are: (x1=smoking, x2=obese) or
#no.yes=c("No","Yes");
#Alternatively, y can be proportion of successes (alive).


#R7.1.6  Chauvinistic attitude to women's work example
library(HSAUR)
data(womensrole); attach(womensrole)
wm=glm(cbind(agree, disagree) ~ sex + education,
       data = womensrole, family = binomial())
su=summary(wm); su
ane=rep(1, length(sex))
bigx=cbind(ane, sex, education)#construct the big X matrix
xbeta=bigx %*% wm$coef #compute X times beta
W=1/(1+exp(xbeta)) #this is the W^ in (7.1.15)
partial=mean(W)*mean(1-W)*wm$coef[2:3]
print("partials of agreement probability w.r.t. sex and
education"); partial
wm2=glm(cbind(agree, disagree) ~ sex +
          education+sex*education,
        data = womensrole, family = binomial())
summary(wm2) #note that education and interaction are
#significant



#R7.1.7 Study of satisfaction with housing High, medium
#or low.
library(MASS)
house.plr <- polr(Sat ~ Infl + Type + Cont, weights=
                    Freq, data = housing);house.plr; summary(house.plr)



#R7.1.8  study of consumer choice between four yogurt brands
#choice= one of yoplait, dannon, hiland, weight (wt watcher)
library(Ecdat); library(VGAM) #Lee's vector glm and gam
data(Yogurt); attach(Yogurt)
nam=names(Yogurt)
nam2=nam[2:8]
Yog2=Yogurt[,nam2] #create data for all variables except id
#in first column and choice in the last column of data
summary(Yog2)  #avoid having to type 7 names exactly
fit = vglm(choice ~ ., multinomial, data=Yog2) #the dot picks
# all choices without having to write them out
coef(fit, matrix=TRUE)
#output matrix lists effect on odds ratio of each regressor



#R7.2.2  Tobit Model Illustration from package "survival"
library(survival)#this package comes with Tobin data
data(tobin) #loads Tobin data for durable expenditure, age,
#and quantity of liquid assets.
#help(tobin)
attach(tobin)
myc=c("D","A","L")
matplot(cbind(durable, age/10,quant/100),pch=myc,typ="b",
        main="Rescaled Tobin data on Durable Goods Purchases",
        xlab="Observation Number", ylab="Durable exp., Age/10,
 Liquidity/100")
tfit <- survreg (Surv(durable, durable>0, type="left")
                 ~age + quant, data=tobin, dist="gaussian")
#(Intercept)         age       quant
#15.14486636 -0.12905928 -0.04554166
#Scale= 5.57254
#type means type of censoring, for Tobin it is left censoring
prs=predict(tfit,type="response")
#pre=tfit$co[1]+tfit$co[2]*age+tfit$co[3]*quant
mysu=as.numeric(Surv(durable, durable > 0, type = "left"))
matplot(cbind(mysu, prs),pch=c("D","F"),typ="b",
        main="Compare Fitted and Actual Durable Goods Purchases",
        xlab="Observation Number",ylab="D=Durable exp.,F=Fitted Exp")
reg1=lm(durable~age+quant)
f2=fitted(reg1)
matplot(cbind(durable, prs,f2),pch=c("D","C","O"),typ="b",
        main="Fitted (Censored or OLS) and Actual Durable Purchases",
        xlab="Observation Number", ylab="D=Durable exp.,C= censored,
 O=OLS")



#R7.3.1  Female Labor force participation using Inverse Mills
#Ratio
library(micEcon)
data( Mroz87 )
#create a new variable called kids from two variables
Mroz87$kids  <- ( Mroz87$kids5 + Mroz87$kids618 > 0 )
greene <- heckit( lfp ~ age + I( age^2 ) + faminc + kids +
                    educ, wage ~ exper + I( exper^2 ) + educ + city, Mroz87 )
summary( greene$probit ) # summary of the 1st step probit
#                           estimation
summary( greene )        # print entire summary
greene$sigma             # estimated sigma
greene$rho               # estimated rho
summary(greene$invMillsRatio) #information on 2-nd regressor


lin="https://faculty.fordham.edu/vinod/oildata.txt"
oil=read.table(file=lin,header=T)
print("Oil and gas CEO survival data")
head(oil,3);tail(oil,3)



#R7.4.1 Plots Weibull-type hazards and fits parametric CEO
#models
layout(matrix(1:3, nrow = 3, ncol = 1)) #will plot 3 figures
plot(dweibull(1:5,shape=0.2, scale=1), typ="l", main="Weibull
Density",ylab="shape=0.02")
plot(dweibull(1:5,shape=0.7, scale=1), typ="l",
     ylab="shape=0.7")
plot(dweibull(1:5,shape=0.99, scale=1), typ="l",
     ylab="shape=0.99")
layout(matrix(1:1, nrow = 1, ncol = 1)) #remove the layout
lin="https://faculty.fordham.edu/vinod/oildata.txt"
oil=read.table(file=lin,header=T)
attach(oil) #make column names accessible variables.
library(rms) #use Harrell's package)
library(survival)
survreg(Surv(tenure, censor) ~ outside+sales+age)
f= psm(Surv(tenure, censor) ~ outside+sales+age,
       dist="weibull", y=TRUE)
f #prints results to the screen
f= psm(Surv(tenure, censor) ~ outside+sales+age,
       dist="loglogistic", y=TRUE)
f  #prints results



#R7.4.2 Cox type analysis of Oil and Gas Company CEO survial
lin="https://faculty.fordham.edu/vinod/oildata.txt"
oil=read.table(file=lin,header=T)
attach(oil); library(survival)
survreg(Surv(tenure, censor) ~ outside+sales+age)
fit=coxph(Surv(tenure, censor) ~ outside+sales+age)
plot(survfit(fit), xlab="Survival years",
     main="Plot of Fitted Survival Curves for Oil&Gas CEOs",
     ylab="Survival probability with confidence limits")


#R7.4.3 CEO Survival Analysis with sophisticated rms
#Package
lin="https://faculty.fordham.edu/vinod/oildata.txt"
oil=read.table(file=lin,header=T)
attach(oil); library(rms) #bring Harrell's package into
#the memory of R
fit2=cph(Surv(tenure,censor)~outside+sales+age,data=oil)
fit2 #print results

#cph is Cox Prop Hazard Model in this package
#Buckley-James distribution-free least squares regression
fit3=bj(Surv(tenure, censor)~outside+sales+age,data=oil)
fit3  #print B-James results,
fitp=psm(Surv(tenure, censor)~outside+sales+age,data=oil)
fitp#print parametric survival model



#R7.4.Out1 OUTPUT of Cox Proportional Hazards Model
# obtained by using the "cph" function of rms library
#Model Tests    Discrimination
#Indexes
#Obs         80    LR chi2     23.29    R2       0.253
#Events      72    d.f.            3    Dxy      0.524
#Center -3.0625    Pr(> chi2) 0.0000    g        0.718
#Score chi2  23.77    gr       2.050
#Pr(> chi2) 0.0000
#coef  se(coef)    z     p-val
#outside  3.86e-02 6.79e-02  0.569 0.569462
#sales    7.18e-05 2.94e-05  2.445 0.014505
#age     -6.14e-02 1.77e-02 -3.478 0.000505



#R7.4.Out2 OUTPUT Buckley-James Regression Model
#by using the "bj" function
#Obs     Events       d.f. error d.f.      sigma
#80         72          3         68     0.7125
#Value Std. Error       Z     Pr(>|Z|)
#Intercept -7.846e-01  7.215e-01 -1.0875 2.768e-01
#outside   -3.612e-02  5.586e-02 -0.6466 5.179e-01
#sales     -3.981e-05  2.033e-05 -1.9584 5.019e-02
#age        5.820e-02  1.063e-02  5.4733 4.417e-08



#R7.4.4 CEO Survival with categorical classifications
# the snippet #R7.4.3 must be in memory for following to work
#define factor (categorical variable) for
#outside directors below the data median of 5 such directors.
fac.outside=factor(outside <5)#define for number of outside
#directors<5
#anova statistically significant effect
survdiff(Surv(tenure,censor)~fac.outside)
fito=survfit(Surv(tenure, censor) ~ fac.outside)
survplot(fito, xlab="Years as the CEO")
rhs=cbind(outside, sales, age)
hazard.ratio.plot(cbind(rhs), Surv(tenure, censor),
                  which=1,cex=1, main="Hazard for CEOs from Outside Directors")



#R7.4.Out3  OUTPUT Chisq test for outside directors
survdiff(formula = Surv(tenure, censor) ~ fac.outside)
#N Observed Expected (O-E)^2/E (O-E)^2/V
#fac.outside=FALSE 48       43     29.9      5.79      12.1
#fac.outside=TRUE  32       29     42.1      4.10      12.1
#
#Chisq= 12.1  on 1 degrees of freedom, p= 0.000501
#Log Hazard Ratio evaluation times: 4.575 19.750
#$log.hazard.ratio
#[,1]          [,2]
#outside  9.036506e-02 -0.0731265607
#sales    7.892914e-05  0.0002626575
#age     -1.344592e-01 -0.0357473767