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


p=seq(from=0.01,to=0.99, by=0.01);p
q=1-p
H=p*log(1/p)+q*log(1/q)
plot(p,H,typ="l")





#R12.8.1a panel data models: Get data
library(plm);data(Grunfeld);attach(Grunfeld)
head(Grunfeld,2);tail(Grunfeld,2)





#R12.8.1b grahics for panel data models
coplot(inv ~ year|firm, type="l")
coplot(inv ~ year|firm, type="b")
car::scatterplot(inv~year|firm, boxplots=FALSE, 
                 smooth=TRUE, data=Grunfeld)
gplots::plotmeans(inv ~ firm, main=
                    "Heterogeineity across firms", data=Grunfeld)
gplots::plotmeans(inv ~ year, main=
                    "Heterogeineity across years", data=Grunfeld)






#R12.8.1c further fixed-effect graphics, panel models
wit=plm(inv~value+capital,data=Grunfeld,model="within")
withat<-fitted(wit)
car::scatterplot(withat~value|firm, boxplots=FALSE,
                 xlab="firm", ylab="fitted inv",smooth=FALSE)
abline(lm(inv~value+capital),lwd=3, col="red")
car::scatterplot(withat~capital|firm,boxplots=FALSE,
                 xlab="firm", ylab="fitted inv",smooth=FALSE)
abline(lm(inv~value+capital),lwd=3, col="red")





#R12.8.1d further panel data models and tests
fixef(wit)#extracts individual fixed effects
twy=plm(inv~value+capital,data=Grunfeld,model="within",
        effect="twoways") #two-way effect time and indivi
fixef(twy, effect="time")#extract time fixed effect
bet=plm(inv~value+capital,data=Grunfeld,model="between")
pool=plm(inv~value+capital,data=Grunfeld,model="pooling")
fdi=plm(inv~value+capital,data=Grunfeld,model="fd")
reff=plm(inv~value+capital,data=Grunfeld,model="random")
rac=pvcm(inv~value+capital,data=Grunfeld,model = "within")
pooltest(wit, pool)# F test fails to implement 
zplm = plm(inv ~ value + capital, data = Grunfeld)
pooltest(zplm, rac)# test does implement






#R12.8.1e  Simpson's paradox avoided by fixed effects
library(plm);library(Ecdat);data(Fatality)
plm.fatality <- pdata.frame(Fatality,index=c("state","year"))
mp=plm(mrall ~ beertax,data=plm.fatality,model="pooling")
coef(mp) #positive slope 0.36
mw=plm(mrall ~ beertax,data=plm.fatality,model="within")
coef(mw) #negative -0.66





#R12.8.2 panel data GMM for Arellano and Bond case
gm=pgmm(dynformula(emp ~ wage + capital +output, lag=list(2,
                                                          1, 0, 1), log = TRUE), EmplUK, effect = "twoways", model = 
          "twosteps", gmm.inst = ~log(emp), lag.gmm = list(c(2, 99)))
summary(gm)