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



#R6.1.1  Two equation Haavelmo model from Klein I data
#load the library, data, and construct Y=C+I using (6.1.1)
#definition
library(sem)
data(Klein)
attach(Klein)#now names C and I are accessible
summary(C); length(C) #T=22 here
Y=C+I  #define the Y data from the definitional Eq. (6.1.1)
reg1=lm(C~Y)  #this is Eq. (6.1.2)  b=0.73968   t-statistic=
#12.429 a= 12.90282    and its t statistic=   3.911


#R6.1.2  Indirect Least squares for the Haavelmo model
# Be sure to run #R6.1.1 in R memory before following this
reg2=lm(C~I); summary(reg2)#intercept=52.1804,slope=0.8782
reg3=lm(Y~I); summary(reg3)#intercept=52.1804,slope=1.878
b=reg2$coef[2]/reg3$coef[2];b  #0.4675748=ILS slope or MPC
a=reg2$coef[1]*(1-b) #27.78216 ILS estimate of intercept



#R6.2.1  Instrumental Variables estimation of consumption
#function. Load the library and data, and construct Y
#by adding consumption and investment.
#We use total production of private industry as W
library(sem)
data(Klein)
attach(Klein)#Now names C and I are accessible.
#X=production of private industry, TPPI, in sem package notation
summary(C)
bigt=length(C)
ane=rep(1,bigt) #T=22 here
bigx=cbind(ane, C+I)  #define the income from the definitional
#        Eq. (6.1.2)
y=C #notation for dependent variable is y in (6.2.1)
bigw=cbind(ane,X) #define the instrument matrix
inver = solve(t(bigw) %*% bigx)
biv= inver %*% (t(bigw) %*% y); biv   #this is the IV estimator
resid=y - bigx %*% biv
sig2 = sum(resid^2)/bigt
vbiv= inver %*% crossprod(bigw) %*% inver  #see Eq. (6.2.7)
se=sqrt(diag(vbiv));se #standard errors
tstat=biv/se; tstat #t stat for IV estimator
#Next, verify that GIVE and IV estimators are identical here
hatw=bigw %*% solve(crossprod(bigw)) %*% t(bigw)
bgive1=solve(t(bigx) %*% hatw %*% bigx)
bgive2=t(bigx) %*% hatw %*% y
bgive=bgive1 %*% bgive2
cbind(biv, bgive)  #identical coefficients


#R6.2.2   Wu\(-\)Hausman test example showing need for IV
#estimator.
#be sure to Run this only after #R6.2.1 is run in the memory
#of R.
reg2=lm(y~bigx+resid-1); summary(reg2) #OLS is inconsistent
#since t stat on resid as regressor is highly significant.
bigy=C+I  #define the offending regressor potentially related
#to errors
resid2=resid(lm(bigy~X))  #regress offending variable on the
#instrument
reg3=lm(y~bigx+resid2-1); summary(reg3) #OLS is inconsistent
#since t stat= -3.042 on resid2 as regressor is still highly
#significant.


#R6.2.3  Data Download and use packages
lin="https://faculty.fordham.edu/vinod/ch6data.txt"
lin <- read_table("C:/Users/reyna/Desktop/libros y sus pdf/Hands-on-intermediate-econometrics-using-r/los txt/ch6data.txt")
ch6data=read.table(file=lin,header=T);head(ch6data,3)
library(generalCorr); options(np.messages=FALSE)
c1=causeSummary(ch6data)
xtable::xtable(c1)



#R6.3.1   2SLS estimator for the Haavelmo model.
library(sem); data(Klein); attach(Klein); Y=C+I
reg2=lm(C~I); summary(reg2)#this is reduced form OLS for C
reg3=lm(Y~I); summary(reg3)#this is reduced form OLS for Y
# next compute fitted values from reduced form estimation
Chat=fitted(reg2)# fitted C is not used
Yhat=fitted(reg3)
#these replace the income on RHS of consum function
reg4=lm(C~Yhat)  #this is second stage OLS
summary(reg4) #these are 2SLS estimators



#R6.3.2 Hands-On Algebra for the k-class estimator.
library(sem); data(Klein); attach(Klein); T=length(C);Yj=C+I
reg3=lm(Yj~I); Yhat=fitted(reg3); Xj=rep(1,T);
Zj=cbind(Yj,Xj)
Zhatj=cbind(Yhat,Xj)
#Yj, Y^j Vj and V^j are all 22 by 1 vectors here
twoby2=t(Zhatj) %*% Zhatj ; twoby2
#         Yhat     Xj
#Yhat  66679.26   1203
#Xj    1203.00     22
Vhatj=Yj-Yhat
#verify orthogonality of fitted values with residuals
t(Vhatj) %*% Yhat # the answer is smaller than 0.00000000000
A=crossprod(Vhatj)
yy=crossprod(Yj)
yy-A  #this equals  66679.26 top left corner of 2x2 matrix
#above.


#R6.3.3 Verify that k-class estimator=2SLS when k=1.
#Run R6.3.2 into memory of R before the following
k=1; mtx11=yy-k*A  #top left corner on the needed matrix
mtx12=t(Yj) %*% Xj
mtx21=t(Xj) %*% Yj
mtx22=t(Xj) %*% Xj
vec11=t(Yj)-k* t(Vhatj)
vec21=t(Xj)
mtx1=cbind(mtx11,mtx12)#create first row of matrix
mtx2=cbind(mtx21,mtx22)#create second row of matrix
mtx=rbind(mtx1,mtx2)
vec=rbind(vec11,vec21)
bkclass=solve(mtx) %*% vec %*% C; bkclass  #implement (6.3.3)
#the intercept =27.782 and b=MPC= 0.4676 agrees with theory



#R6.3.4 LIML estimator as k-class estimator for the Haavelmo
#model
#Run R6.3.2 and R6.3.3 into memory of R before the following
hat0= bigx %*% solve(crossprod(bigx)) %*% t(bigx)
hatj= Xj %*% solve(crossprod(Xj)) %*% t(Xj)
#hatj is 22 x 22 degenerate matrix with all
#elements=0.04545455
Y=C+I; Y0j=cbind(C,Y)
E0j=(diag(T)-hat0) %*% Y0j
E1j=(diag(T)-hatj) %*% Y0j
W0j=t(E0j) %*% E0j
W1j=t(E1j) %*% E1j
W=solve(W1j) %*% W0j
eiw=eigen(W); eiw$values
lamda.min=min(eiw$values)  #lamda min is zero due to
#degeneracy
E0jb=C-mean(C)
E1jb=resid(lm(C~I))
W0jb=t(E0jb) %*% E0jb
W1jb=t(E1jb) %*% E1jb
ev=W0jb / W1jb; ev  #1.209153 is a sensible value for
#k-class k
mtx11=yy-ev*A
vec11=t(Yj)-ev* t(Vhatj)
vec21=t(Xj)
mtx1=cbind(mtx11,mtx12)#create first row of matrix (6.3.3)
mtx2=cbind(mtx21,mtx22)#create second row of matrix
mtx=rbind(mtx1,mtx2)
vec=rbind(vec11,vec21)
bLIML=solve(mtx) %*% vec %*% C; bLIML
#slope=0.5984009  Intercept= 20.6283499



#R 6.4.1 Kmenta demand/supply estimation by various models
rm(list=ls()) #rm means remove them all to clean up old stuff
library(systemfit); data("Kmenta"); attach(Kmenta)
eqDemand = consump ~ price + income #demand eq exo income
eqSupply = consump ~ price + farmPrice + trend #supp eq
eqSystem = list(dd = eqDemand, sup = eqSupply)
#above defines the equation system with a list
fitols = systemfit(eqSystem)
#when no method is specified, it gives OLS
print(fitols)
b.ols=coef(fitols)
fit2sls = systemfit(eqSystem, method = "2SLS", inst = ~income +
                      farmPrice + trend, data = Kmenta)
b.2sls=coef(fit2sls);summary(fit2sls)
restrict = "dd_price + sup_farmPrice = 0"
fit2slsRmat = systemfit(eqSystem, method = "2SLS",
                        inst = ~income + farmPrice + trend, data = Kmenta,
                        restrict.matrix = restrict)
b.restr.2sls=coef(fit2slsRmat);summary(fit2slsRmat)
options(digits=4) #allows compact printing
rbind(b.ols, b.2sls,b.restr.2sls)
UnRSS=sum(resid(fit2sls)^2)# 162.4
ReRSS=sum(resid(fit2slsRmat)^2)# 161.7
df=summary(fit2sls)$df[2]  #7, 33
#summary(fit2slsRmat)$df  # 6, 34

r=1
#Fstat=((ReRSS-UnRSS)/r ) /(UnRSS/df );Fstat
lrtest(fit2sls, fit2slsRmat)



#R6.5.1 Wold Recursive Model Cobb-web agricultural supply model
library(systemfit)
data( "Kmenta" ); attach(data.frame(Kmenta)); q=consump
#change of notation q for supply
reg1=lm(q ~ farmPrice + trend); summary(reg1)
reg2=lm(price ~q+income); summary(reg2)



#R6.5.2 Demand Supply model illustrates Kronecker product &
#inverse
# Run all commands of #R6.4.1 till fit2sls in memory of R
a=matrix(c(1,3),1,2);b=matrix(10:13,2,2);a;b; a %x% b
#above a is 1x2 and b is 2x2 and Kronecker product is 2x4
#  10   12   30   36
#  11   13   33   39
omega.hat =cov(resid(fit2sls));omega.hat
#cov of 20 by 2 matrix
#         demand   supply  This is omega.hat matrix
#demand 3.459426 3.782355
#supply 3.782355 5.085960
I20=diag(20) # create identity matrix of size 20 in R
eed=omega.hat %x% I20  #this is a huge 40 x 40 matrix, not
#printed
max(solve(eed)- solve(omega.hat) %x% I20)  #this is zero


#R6.5.3 Demand Supply model test showing 3SLS is efficient
library(systemfit); data("Kmenta"); attach(Kmenta)
eqDemand = consump ~ price + income #demand eq exo income
eqSupply = consump ~ price + farmPrice + trend
eqSystem = list(dd = eqDemand, sup = eqSupply)
fit2sls = systemfit(eqSystem, method = "2SLS", inst = ~income +
                      farmPrice + trend, data = Kmenta)
fit3sls = systemfit(eqSystem, method = "3SLS",
                    inst = ~income + farmPrice + trend)
b.3sls=coef(fit3sls);summary(fit3sls)
fitI3sls = systemfit(eqSystem, method = "3SLS",
                     inst = ~income + farmPrice + trend, data = Kmenta,
                     maxit = 250)
b.Iter3sls=coef(fitI3sls);summary(fitI3sls)
## perform the Hausman test
h <- hausman.systemfit( fit2sls, fit3sls )
print( h )#Hausman = 2.5357, df = 7, p-value = 0.9244
#p-value being large we accept the null, 3SLS is good!



#R6.5.4 Klein I model with systemfit package.
library(systemfit)#type
#help(package=systemfit) for latest info.
data( "KleinI" )
eqConsump  = consump ~ corpProf + corpProfLag + wages
eqInvest   = invest ~ corpProf + corpProfLag + capitalLag
eqPrivWage = privWage ~ gnp + gnpLag + trend
inst = ~ govExp + taxes + govWage + trend + capitalLag +
  corpProfLag + gnpLag
system = list( Consumption = eqConsump, Investment = eqInvest,
               "Private Wages" = eqPrivWage )
# OLS
kleinOls = systemfit( "OLS", system, data = KleinI )
summary( kleinOls )
# 2SLS
klein2sls = systemfit( "2SLS", system, inst = inst,
                       data = KleinI, rcovformula = 0 )
summary( klein2sls )
# 3SLS
klein3sls = systemfit( "3SLS", system, inst = inst,
                       data = KleinI, rcovformula = 0 )
summary( klein3sls )
# I3SLS
kleinI3sls = systemfit( "3SLS", system, inst = inst,
                        data = KleinI, rcovformula = 0, maxit = 500 )
summary( kleinI3sls )