# R 3.0.1
## Ex5_panel_data_R.txt
## panel data ONLY for yearly/consecutively indexed data
  
library("AER")
library("plm")
data(Grunfeld, package = "AER")
## firm        ... name of firm
## inv(invest) ... investment expenditures
## value       ... value of firm
## capital     ... stock of plant and equipment
Grunfeld                                             # list the data
## writing/reading of data
#write.table(Grunfeld, file = "Grunfeld.txt", col.names = TRUE) 
#dat <- read.table("Grunfeld.txt", header = TRUE)    # read data from txt

gr <- Grunfeld

## choose a subset of the data
#gr <- subset(gr, firm %in% c("General Electric", "General Motors", "IBM") )
##(1)
pgr <- plm.data(gr, index = c("firm", "year")) # define panel-data-structure 
   ## data are in stacked form: [(firm1,t=1,...,T) (firm2, t=1,...,T) ...]'
pgr.n <- 11;  pgr.T <- 20;            # 11 firms, 1935-1954

##(2)
##pooled regression (overall intercept automatically included)
## without overall intercept: ... invest ~ value + capital + 0, ...
gr_pool <- plm(invest ~ value + capital, data = pgr, model = "pooling")
summary(gr_pool)

##(3)
##fixed effects panel data estimation:  within,  individual intercepts
#R2 of plm-within = R2 'without indiv intercepts'
gr_fe <- plm(invest ~ value + capital, data = pgr, model = "within", effect = "individual")
summary(gr_fe)
summary( fixef(gr_fe, effect= "individual") )  #individual intercepts OR fixef(gr_fe)

##(4)
# F-test for the group effects wrt to the simple pooling with one intercept
##(4a)
r2_pool <- summary(gr_pool)$r.squared[1]
#OR:  r2_pool <- cor(pgr$invest, pgr$invest - gr_pool$residuals)^2
##(4b)   #gr_fe$residuals are the fe-residuals
r2_fe <- cor(pgr$invest, pgr$invest - gr_fe$residuals)^2   
df1 <- pgr.n - 1;     df2 <- pgr.n*pgr.T - pgr.n - 2; 
f_fe_pool <- ( (r2_fe - r2_pool) / df1 ) / ( (1 - r2_fe) / df2 )
p_value <- 1- pf(f_fe_pool, df1, df2)
print( c(df1,df2,as.numeric(f_fe_pool), as.numeric(p_value) ) )     
 # OR by    H0: pool    HA: FE
pFtest(gr_fe, gr_pool)                             

##(5)
##fixed effects panel data estimation:  within,  time intercepts 
gr_fe_t <- plm(invest ~ value + capital, data = pgr, model = "within", effect = "time")
summary(gr_fe_t)
summary( fixef(gr_fe_t, effect= "time") ) # time intercepts  OR  fixef(gr_fe_t)   

#between (omitted)
#gr_fe_bi <- plm(invest ~ value + capital, data = pgr, model = "between", effect = "individual")
#gr_fe_bt <- plm(invest ~ value + capital, data = pgr, model = "between", effect = "time")
# fixef only for within 

##(6)
##fixed effects panel data estimation:  within,  both individual and time intercepts 
gr_fe2 <- plm(invest ~ value + capital, data = pgr, model = "within", effect = "twoways")
summary(gr_fe2)
summary( fixef(gr_fe2 , effect = "individual") ) 
summary( fixef(gr_fe2 , effect = "time") )      # NO fixef(gr_fe2 , effect = "twoways")  

# alternatively  Lagrange multiplier test: pool vs fe (indiv, time eff or twoways)
# (omitted)
#plmtest(gr_pool, effect = "individual")
#plmtest(gr_pool, effect = "time")                 # here hardly siginificant
#plmtest(gr_pool, effect = "twoways")

##(7)
##random effects
gr_re <- plm(invest ~ value + capital, data = pgr, model = "random", random.method = "walhus")
summary(gr_re)                # alpha_i ... individual, u_{i,t} ... idiosyncratic

##(8)
## Hausman test:   H0: RE vs HA: FE 
phtest(gr_fe, gr_re)

##(9)
# Durbin-Watson-Test for panel data
pdwtest(gr_fe)

###################
## dynamic panel ##
###################
## data requirement: ONLY small T, T < N do work ##
##(10)
data(Grunfeld, package = "plm")
ppgr <- Grunfeld 
ppgr1 <- subset(ppgr, year %in% c(1935:1940) )

# formula/model
form <- inv ~ value + capital

# Arellano-Bond GMM
gr.ab <- pgmm( dynformula(form, list( 1, c(0,1), 0 )), 
   data = ppgr1, 
   effect = "individual", model = "onestep", 
   gmm.inst = ~ inv , lag.gmm = list(c(2, 3)) )
summary(gr.ab)

# OR identically, in a different notation
gr.ab1 <- pgmm( inv ~ lag(inv, 1) + value + lag(value, 1) 
           + capital | lag(inv, 2:3 ) ,
            data = ppgr1, index = c("firm", "year"),  
            effect = "individual", model = "onestep")
summary(gr.ab1)

##(11)
## dynamic panel EmplUK
data("EmplUK", package ="plm") ## firm: 1-140,  year: 1977-1983
# panel-data-structure already defined
#pem <- plm.data(EmplUK, index = c("firm", "year"))
  
# formula / model
form <- log(emp) ~ log(wage) + log(capital) + log(output)

# Arellano-Bond GMM      effect = "twoways"
empl.ab <- pgmm( dynformula(form, list(2, 1, 0, 2)), 
   data = EmplUK, index = c("firm", "year"), 
   effect = "twoways", model = "twosteps", 
   gmm.inst = ~ log(emp), lag.gmm = list(c(2, 99)) )
summary(empl.ab)

# Arellano-Bond GMM     effect = "individual",
empl1.ab <- pgmm( dynformula(form, list(1, 1, 0, 0)), 
   data = EmplUK, index = c("firm", "year"), 
   effect = "individual", model = "twosteps", 
   gmm.inst = ~ log(emp), lag.gmm = list(c(2, 99)) )
summary(empl1.ab)
#OR identically,  
# Arellano-Bond GMM
empl2.ab <- pgmm( dynformula(form, list(1, c(0,1), 0, 0)), 
   data = EmplUK, index = c("firm", "year"), 
   effect = "individual", model = "twosteps", 
   gmm.inst = ~ log(emp), lag.gmm = list(c(2, 99)) )
summary(empl2.ab)