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



#R9.1.1 Cobb Douglas regression for Metals data already in R,
library(Ecdat); data(Metal) #AS in snippet  #R1.2.1
met=as.matrix(Metal); y=met[,1]; L=met[,2]; K=met[,3]
Ly=log(met[,1]); LL=log(met[,2]); LK=log(met[,3])
reg1 = lm(Ly~LK+LL); summary(reg1)
a1=reg1$coef[1]
b1=reg1$coef[2]
c1=reg1$coef[3]
##  Now nonlinear least squares
reg=nls(formula=(y~a*(L^b)*(K^c)), start=list(a=a1,b=b1,c=c1))
su=summary(reg); su
a=su$coef[1]
b=su$coef[2]
c=su$coef[3]
se.a=su$coef[1,2]
se.b=su$coef[2,2]
se.c=su$coef[3,2]
a #estimated alpha =  2.736079
se.a #Std. Error of est of alpha =  0.92121
b #estimated beta =  0.4036325
se.b #Std. Error of beta =0.1666741
c #estimated gamma = 0.5508710
se.c #Std. Error of c=0.1331163
SCE=b+c  #scale elasticity from nonlinear estimation
v=vcov(reg)#covariance matrix of coefficients, nonlinear model
se=sqrt(se.b^2+se.c^2+2*v[2,3])#SE of SCE using covariance term
v[2,3] # cov(beta,gamma)= -0.0212709 is negative
SCE #estimated scale elasticity = 0.9545035
se #Std. Error of scale elasticity =  0.05439119
#asymptotic normal 95% confidence interval
taul=SCE-1.96*se
taul #lower bound is 0.8478968
tauu=SCE+1.96*se
tauu #upper bound is # 1.061110


#R9.2.1 parametric bootstrap for Cobb\(-\)Douglas Metals model

library(Ecdat); data(Metal) #AS in snippet  #R1.2.1
met=as.matrix(Metal); y=met[,1]; L=met[,2]; K=met[,3]
Ly=log(met[,1]); LL=log(met[,2]); LK=log(met[,3])
reg1 = lm(Ly~LK+LL); summary(reg1)
bigt=length(Ly)
df=summary(reg1)$df[2] #degrees of freedom extracted
ane=rep(1,bigt)#create column of ones
bigx=cbind(ane,LK,LL) #X matrix
htt= hat(bigx) #hat matrix diagonals
bigj=999
se=(summary(reg1)$sigma)*sqrt(ane-htt)
set.seed(2345)
thet=rep(NA,bigj)# place to store theta
par(mfrow=c(2,1)) # two histograms in one plot
for (j in 1:bigj) {
  resj=rnorm(bigt)*se  #Step 1
  yj=bigx%*%reg1$coef + resj  #Step 2
  bj=lm(yj~bigx-1) #Eq. (9.2.10) implemented here
  thet[j]=bj$coef[2]+bj$coef[3]
}  #end of simulation loop for j

hist(thet, main="Histogram for simulated Normal SCE, Metals
data",
     xlab="scale elasticity adding two coefficients")
sort.thet=sort(thet)
LOW=sort.thet[25]# 0.8645969  is lower limit by (9.2.11)
UPP=sort.thet[975]# 1.088195 is upper limit by (9.2.11)
LOW; UPP;
ols.thet=reg1$coef[2]+reg1$coef[3]
LOW2=2*ols.thet -UPP; LOW2#0.8692241 lower limit by (9.2.12)
UPP2=2*ols.thet -LOW; UPP2#1.092822 upper limit by (9.2.12)
#Now simulate with student's t
set.seed(237)
thet=rep(NA,bigj)# place to store theta
for (j in 1:bigj) {
  resj=rt(n=bigt, df=df)*se  #Step 1
  yj=bigx%*%reg1$coef + resj  #Step 2
  bj=lm(yj~bigx-1) #Eq. (9.2.10) implemented here
  thet[j]=bj$coef[2]+bj$coef[3]
}  #end of simulation loop for j
hist(thet, main="Histogram for simulated t-density SCE,
Metals data",
     xlab="scale elasticity adding two coefficients")
sort.thet=sort(thet)
LOW=sort.thet[25]#   is lower limit by (9.2.11)
UPP=sort.thet[975]#  is upper limit by (9.2.11)
LOW; UPP; #0.8601696 and  1.102989
ols.thet=reg1$coef[2]+reg1$coef[3]
LOW2=2*ols.thet -UPP; LOW2#0.8544301 Lower limit by (9.2.12)
UPP2=2*ols.thet -LOW; UPP2#1.097249 Upper limit by (9.2.12)
par(mfrow=c(1,1))
plot(density(thet), main="Histogram for simulated t-density
SCE, Metals data", xlab="scale elasticity adding two
coefficients")



#R9.3.1 Plotting ECDF for Metals data based on residuals
# You must run snippet #R9.1.1 before the following
par(mfrow=c(2,1)) #begin plots of ECDF for residuals
res1=resid(reg1)# recall the residuals of lm model
ecdres=ecdf(res1)# R function ecdf does all for ECDF
plot(ecdres, verticals= TRUE, do.p = FALSE, main="Metals data
ECDF of residuals", xlab="Log-linear model residual",
     ylab="Cumulative probability")
res=resid(reg);ecdre=ecdf(res)
plot(ecdre, verticals= TRUE, do.p = FALSE, main="Metals data
ECDF of residuals", xlab="Nonlinear model residual",
     ylab="Cumulative probability")
mres=mean(res)#9.168556
sdres=sd(res)#448.8478
set.seed(340); xj=rnorm(27, mean=mres, sd=sdres)
ecdxj=ecdf(xj)#empirical CDF values from xj
par(mfrow=c(2,1)) #begin first of two plots
plot(ecdxj, verticals= TRUE, do.p = FALSE, main="Metals data
ECDF of parametric normal residuals", xlab="Nonlinear model
random residual", ylab="Cumulative probability")
xtj=rt(n=27, df=24, ncp=-mres/sdres)
ecdxtj=ecdf(xtj)
plot(ecdxj, verticals= TRUE, do.p = FALSE, main="Metals data
ECDF of parametric t residuals", xlab="Nonlinear model random
residual", ylab="Cumulative probability")



#R9.3.2 non-parametric iid bootstrap for Metals using residuals
library(Ecdat); data(Metal) #AS in snippet  #R1.2.1
met=as.matrix(Metal); y=met[,1]; L=met[,2]; K=met[,3]
Ly=log(met[,1]); LL=log(met[,2]); LK=log(met[,3])
reg1 = lm(Ly~LK+LL); summary(reg1)
library(simpleboot); set.seed(345)
bootiid=lm.boot(reg1, R=999)
summary(bootiid)  #prints bootstrap standard errors for coef.s
#plot(bootiid)  works only for 1 regressors
s <- samples(bootiid, "coef")
## Histogram for the intercept
hist(s[1,])
perc.lm(bootiid, p=c(0.025,0.975))  #gives a 95 % conf. int.



#R9.3.3  Metals data nonlinear model SINGLE BOOTSTRAP
rm(list=ls()); date()  #clean the memory of R
library(Ecdat); data(Metal) #AS in snippet  #R1.2.1
met=as.matrix(Metal); y=met[,1]; L=met[,2]; K=met[,3]
Ly=log(met[,1]); LL=log(met[,2]); LK=log(met[,3])
reg1 = lm(Ly~LK+LL); summary(reg1)
a1=reg1$coef[1]; b1=reg1$coef[2]
c1=reg1$coef[3]
reg=nls(formula=(y~a*(L^b)*(K^c)), start=list(a=a1,b=b1,c=c1))
su=summary(reg)
a=su$coef[1];b=su$coef[2];c=su$coef[3]
SCE=b+c  #This is just scale elasticity
se.a=su$coef[1,2];se.b=su$coef[2,2];se.c=su$coef[3,2]
v=vcov(reg)#covariance matrix of coefficients, nonlinear model
se=sqrt(se.b^2+se.c^2+2*v[2,3])#SE of SCE
B=999 #size of single bootstrap
n=length(su$residuals) #number of observations
# Notation in the code is PTS=place to store
bintercept2=rep(0,B)#PTS estimated alpha from stage 1 boot
bL2=rep(0,B) #PTS estimated beta from each  stage 1 boot
bK2=rep(0,B) #PTS estimated gamma from each  stage 1 boot
se.b2=rep(0,B) #PTS Std. Error of beta from  stage 1 boot
se.c2=rep(0,B) #PTS SE(gamma) from each  stage 1 boot
cov.bc2=rep(0,B) #PTS cov(beta,gamma) from  stage 1 boot
errors=matrix(0,n,B) #PTS residuals from each  stage 1 boot
resid=su$residuals #residuals from initial nonlinear regr.
# Begin boostrap resampling with random seed  ###  ### ##
set.seed(239) #sample(.) uses the seed
for (j in 1:B){
  res2=sample(x=resid, size=n, replace=T)
  QSim2=a*(L^b)*(K^c)+res2
  reg2=nls(formula=(QSim2~a*(L^b)*(K^c)),
           start=list(a=a1,b=b1,c=c1))
  su2=summary(reg2)
  bintercept2[j]=su2$coef[1]#pick stage 1 alpha
  #coeff and save it
  bL2[j]=su2$coef[2]  #pick & save stage 1 beta coef
  bK2[j]=su2$coef[3] #pick  & save stage 1 gamma coef
  se.b2[j]=su2$coeff[2,2] #pick & save the Std. Error of
  #first-stage beta
  se.c2[j]=su2$coeff[3,2]  #pick &save SE(stage1 gamma)
  v2=vcov(reg2)
  cov.bc2[j]=v2[2,3] #pick &save cov(stage1 beta,
  # and first-stage gamma)
  errors[,j]=su2$residuals #pick & save resid stage 1
}  #end the big loop over B (large number like J=999)
##Part 2  Following implements Hall's suggestions
# without using the quantile() function
mean(errors)#
SCE2=bL2+bK2 #999 bootstrap estimates of SCE scale elas.
se2=sqrt(se.b2^2+se.c2^2+2*cov.bc2) # and its standard error
R1=(SCE2-SCE)/se2
n25=round((B+1)*0.025,0)
n975=round((B+1)*0.975,0)
#percentile bootstrap conf int for the single boot
sort.SCE2=sort(SCE2)
sort.SCE2[n25] #lower bound is  0.8503211
sort.SCE2[n975] #upper bound is 1.058657
#following is tempting but Hall has proved it to be wrong
sort.se2=sort(se2)
sort.R1=sort(R1)
SCE-sort.R1[n975]*sort.se2[n975] #lower bound is 0.8128355
SCE-sort.R1[n25]*sort.se2[n25] #upper bound is  1.021643
#percentile-t interval for the single bootstrap
#Hall shows SE from nonlinear estimate is more correct
SCE-sort.R1[n975]*se #lower bound is 0.837974
SCE-sort.R1[n25]*se #upper bound is  1.066788



#R9.3.4 Wild bootstrap using rbinom
set.seed(1415); n=25; x=rnorm(n)   #predictor variable
y=3 + 0.5*x + rnorm(n, sd=2)       #response
My=lm(y~x)                    #linear model fit
MyCoef=coef(My)          #obtain original model coefficients
MyMatrix=model.matrix(My)#obtain the predictor model matrix
Residuals=residuals(My)  #replace with the residuals
##The loop would be set up around this part of the code##
Ans=rep(NA,999)
for (j in 1:999){
  WildRV=rbinom(n, size=1, prob=0.5) #generate a binomial R.V.
  WildRV[WildRV==0]=-1 #Update 0's as -1 as in Rademacher
  WildResi=Residuals*WildRV  #generate wildboot residuals
  WildResponse=MyMatrix %*% as.matrix(MyCoef,nrow=1)+WildResi
  WildBootCoef=coef(lm(WildResponse~x)) #get OLS coef.
  Ans[j]=WildBootCoef[2] }#end j loop
quantile(Ans, probs=c(0.025, 0.975))



#R9.4.1 Implementing Double Bootstrap Metals data including
#GNR double bootstrap. Parts are written by Diana Rudean.
rm(list=ls()); date() #clean up memory at start
library(Ecdat); data(Metal) #AS in snippet  #R1.2.1
met=as.matrix(Metal); y=met[,1]; L=met[,2]; Kap=met[,3]
#RESCALE makes similar units of all X, so nls reliable
y=y/1000; L=L/100; Kap=Kap/1000
Ly=log(y);LL=log(L);LK=log(Kap)
reg1=lm(Ly~LL+LK);
a1=reg1$coef[1]; b1=reg1$coef[2]
c1=reg1$coef[3]
# # NonLinear Regression using nls function with ols
#starting values
reg=nls(formula=(y~a*(L^b)*(Kap^c)),
        start=list(a=a1,b=b1,c=c1))
su=summary(reg)
a=su$coef[1]; b=su$coef[2]
c=su$coef[3]
se.b=su$coef[2,2]; se.c=su$coef[3,2]
SCE=b+c     #   scale elasticity from NLS model
v=vcov(reg)
se=sqrt(se.b^2+se.c^2+2*v[2,3])#  SE(SCE)
resid=su$residuals
n=length(su$residuals)

#J=1999; K=250 choice in McCulough\(-\)Vinod J1998b

J=999  #Bootstrap size denoted as J in Vinod's publications
K=135
n25=round((J+1)*0.025); if (n25<1) n25=1
n975=round((J+1)*0.975)
bintercept2=rep(0,J); bL2=rep(0,J)
bK2=rep(0,J); se.b2=rep(0,J)
se.c2=rep(0,J); cov.bc2=rep(0,J)
errors=matrix(0,n,J); R1=rep(NA, J)
SCE1=rep(NA, J)
bL3=matrix(0,K,J); bK3=matrix(0,K,J)
se.b3=matrix(0,K,J); se.c3=matrix(0,K,J)
cov.bc3=matrix(0,K,J)
set.seed(239)
for (j in 1:J){      #BEGIN  j loop
  res2=sample(x=resid, size=n, replace=T)
  ySim2=a*(L^b)*(Kap^c)+res2
  reg2=nls(formula=(ySim2~a*(L^b)*(Kap^c)),
           start=list(a=a1,b=b1,c=c1))
  su2=summary(reg2)
  bintercept2[j]=su2$coef[1]
  bL2[j]=su2$coef[2]; bK2[j]=su2$coef[3]
  se.b2[j]=su2$coeff[2,2]
  se.c2[j]=su2$coeff[3,2]
  v2=vcov(reg2); cov.bc2[j]=v2[2,3]
  errors[,j]=su2$residuals  #to be resampled in double bootstrap
  SCE1[j]=bL2[j]+bK2[j]
  se1=sqrt(se.b2[j]^2+se.c2[j]^2+2*cov.bc2[j])
  R1[j]=(SCE1[j]-SCE)/se1
}  #end loop for j
mean(errors)
#################################  DOUBLE BOOT
# store all j k results
for (j in 1:J){
  for (k in 1:K){
    res3=sample(x=errors[,j], size=n, replace=T)
    ySim3=bintercept2[j]*(L^(bL2[j]))*(Kap^(bK2[j]))+res3
    
    #first step in McCullough\(-\)Vinod page 87
    #initialize the procedure using nonlinear regr values
    #set GNR estimates a2,b2,c2, y.gnr equal to a,b,c, y
    #from NonLinear regression
    a2=a; b2=b; c2=c
    y.gnr=a2*(L^b2)*(Kap^c2)  #Get pseudo y values for GNR
    #Steps 2 to 4 McCullough\(-\)Vinod page 88
    m=4 #each nonlinear estimation of bootstrap is replaced by
    #4 linear estimations
    
    for (i in 1:m){
      f1=(L^b2)*(Kap^c2) #derivative with respect to alpha (a2)
      f2=a2*LL*(L^b2)*(Kap^c2) #derivative with respect to beta (b2)
      f3=a2*(L^b2)*LK*(Kap^c2) #derivative with respect to gamma (c2)
      r=ySim3-y.gnr  #dependent variable for GNR is residual
      
      #Step 5 McCullough\(-\)Vinod page 88, Auxiliary regression of GNR
      
      reg21=lm(r~f1+f2+f3-1) #regress r on f1, f2, f3
      su21=summary(reg21)
      a21=su21$coef[1]  #coefficients of auxiliary regression
      b21=su21$coef[2]; c21=su21$coef[3]
      
      #update the coefficient estimates  Step 6 McCullough\(-\)Vinod
      
      a2=a2+a21  #sequential updating by GNR coefficients
      b2=b2+b21; c2=c2+c21
      y.gnr=a2*(L^b2)*(Kap^c2) #update the fitted values
    } #END LOOP over i we update 4 times ending GNR updating
    #the standard errors for the coefficients are those from
    #the final regression of r on f1, f2, and f3
    se.b3[k,j]=su21$coef[2,2]
    se.c3[k,j]=su21$coef[3,2]
    v=vcov(reg21); cov.bc3[k,j]=v[2,3]
    bL3[k,j]=b2; bK3[k,j]=c2}  #end k loop
}   #end j loop
##############################  Process Double Boot Results
#bootstrap estimate of SCE returns to scale parameter
#and its standard error (nonlinear bootstrap)

U=rep(NA,J)  #make it NA vector of length J for each k
R2=rep(NA,J); xx=rep(NA,K)
for (j in 1:J){   #LOOP for j
  for (k in 1:K){  #LOOP for k
    SCE2=bL3[k,j]+bK3[k,j]
    se2=sqrt(se.b3[k,j]^2+se.c3[k,j]^2+2*cov.bc3[k,j])
    R2[j]=(SCE2-SCE1[j])/se2
    xx[k]=(R2[j]-R1[j])
  }  #end k loop
  U[j]=length(xx[xx<=0])/K #Beran's variable should be Uniform
}   #end j loop
#percentile bootstrap confidence interval (nonlinear bootstrap)
sort.U=sort(U); sort.U[n25]; sort.U[n975]
low=round(sort.U[n25]*(J+1))
if (low<1) low=1
upp=round(sort.U[n975]*(J+1))
if (upp>J) upp=J
layout(matrix(1:2, nrow = 2, ncol = 1))

hist(U, main="Histogram for U, Gauss\(-\)Newton Regression,
Metals",
     xlab="U = #(R** < or = R*) / K")
##second stage GNR interval
sort.R1=sort(R1)
SCE+sort.R1[low]*se #lower bound is 0.8353043
SCE+sort.R1[upp]*se #upper bound is 1.080800
plot(density(U),main="Density for U, Gauss\(-\)Newton Regression,
Metals",
     
     xlab="U = #(R** < or = R*) / K")
date()  #Helps count the CPU time taken by the bootstrap
#### ############ Testing Uniformity  #### #### ####
#test statistic for uniformity
Uposi=U[abs(U)>0] #consider nonzero U values only
df=2*length(Uposi)
test=(-2)*sum(log(Uposi))#2086.685
test # 1970.315
qchisq(0.95,df=df) #2086.685
pchisq(test,df) #p-value 0.4303943



#R9.5.1  OLS applied to Keynesian consumption function.
library(meboot); library(car); library(lmtest)
data(USconsum); attach(USconsum)
lc <- log(consum); ly <- log(dispinc)#take Logs
lmcf <- lm(lc[2:51] ~ lc[1:50] + ly[1:50])
#cf=consumption function
coeftest(lmcf)
durbin.watson(model = lmcf, max.lag = 4)#higher order DW test



#R9.5.2  meboot applied to consumption function.
#Let entire snippet #R9.5.1 be in R memory before this.
semx = meboot(x = ly, reps = 999, reachbnd = TRUE,
              expand.sd= TRUE,force.clt = TRUE)$ensemble
#Above two lines get 999 replicates of log y series
layout(matrix(1:2, nrow = 2, ncol = 1))
nam=c("original data","j=1","j=2","j=3")
#nam=as.character(1:4)
matplot(semx[,1:4],typ="l", lty=c(1:4), main="US income data
and 3 meboot resamples",  xlab="time 1948 to 1998",
        ylab="log disposable income")
legend(x=2.3,y=3.1,legend=nam, lty=c(1:4), col=c(1:4))
lmrcf <- lm(lc[2:51] ~ lc[1:50])  #restricted OLS of (9.5.5)
bb <- coef(lmrcf)
resi <- resid(lmrcf) #resid restricted model needed below
bigJ <- 999
et <- meboot(x = resi, reps = 1)$ensemble
AC <- arima.sim(n = length(lc), model = list(order = c(1, 0,
                                                       0), ar = bb[2]), n.start = 1, innov = c(lc[1], et + bb[1]),
                start.innov = 0) #next 2 lines gets resamples of y
semy <- meboot(x = AC, reps = bigJ, reachbnd = TRUE,
               expand.cd= TRUE,force.clt = TRUE)$ensemble
matplot(semy[,1:4],typ="l", lty=c(1:4), main="US consumption
data and 3 meboot resamples",  xlab="time 1948 to 1998",
        ylab="log consumption")
legend(x=2,y=3.05,legend=nam, lty=c(1:4), col=c(1:4))



#R9.5.3 "meboot" applied to consumption function.
#Let entire snippet #R9.5.1 and 2 be in the R memory before
#using this.
olsHALL.b <- function(x, y) {
  n <- length(x)
  x <- cbind(y[-n], x[-n])
  y <- y[-1]; tol <- 1e-07
  p <- ncol(x); ny <- NCOL(y)
  lmh <- lm(y ~ x)
  coef(lmh)[3]} #function ends here
bigJ <- 999; bigK <- 100
allK.artif.b <- rep(0, (bigJ * bigK + 1))
all.LO <- rep(0, bigK); all.UP <- all.LO
set.seed(321)
for (k in 1:bigK) {
  et <- meboot(x = resi, reps = 1)$ensemble
  AC <- arima.sim(n = length(lc), model = list(order = c(1,
                                                         0, 0), ar = bb[2]), n.start = 1,
                  innov = c(lc[1], et + bb[1]),
                  start.innov = 0)
  semy <- meboot(x = AC, reps = bigJ,
                 reachbnd = TRUE, expand.sd = TRUE,
                 force.clt = TRUE)$ensemble
  semy <- cbind(AC, semy)
  semx <- meboot(x = ly, reps = bigJ,
                 reachbnd = TRUE, expand.sd = TRUE,
                 force.clt = TRUE)$ensemble
  semx <- cbind(ly, semx)
  all.artif.b <- rep(NA, ncol(semy))
  for (h in seq(along = all.artif.b))
    all.artif.b[h] <- olsHALL.b(x = semx[,h], y = semy[, h])
  if (k == 1)
    allK.artif.b[1] <- all.artif.b[1]
  k1 <- (k - 1) * bigJ + 2
  k2 <- k * bigJ + 1
  allK.artif.b[k1:k2] <- all.artif.b[2:(bigJ + 1)]
  zci <- zero.ci(all.artif.b)
  all.LO[k] <- zci$lolim
  all.UP[k] <- zci$uplim
  if ((k * 100/bigK)%%5 == 0)
    cat(paste((k * 100/bigK), "%",
              sep = ""), "complete.\n")  }
##  this may take a long time
##
## look for the message: 100% complete.
avLO <- mean(all.LO); avLO  #-0.083
avUP <- mean(all.UP); avUP#0.333
# now size corrected sc bounds
cc <- rbind(as.matrix(all.LO), as.matrix(all.UP))
zcc <- zero.ci(cc)  #a function called zero.ci is part of
# the meboot package
scLO <- zcc$lolim; round(scLO,3) #-0.110
scUP <- zcc$uplim; round(scUP,3) #0.396



# R9.5.4 snippet
library(meboot);    set.seed(12345)
xt = rnorm(1:100, mean = 9, sd = 12)
# Vectorize mebootSpear
mebootSpear_vec <- Vectorize(mebootSpear, vectorize.args = "setSpearman")
# Generate 1 replica of the original series xt with the desired rho
Z = do.call(cbind, mebootSpear_vec(xt, reps = 1,
                                   setSpearman = seq(-1,1,.01), drift = FALSE)[2,])
# Expand to 50 standard meboot reps of each rho value
new_MCS = list()
for(i in 1:dim(Z)[2]){new_MCS[[i]] = mebootSpear(Z[,i], reps = 50,
                                                 setSpearman = 1, drift = FALSE)$ensemble}
hist(cor(do.call(cbind,new_MCS), xt),
     main = "mebootSpear Simulation 10,000 Iterations",xlab = "Correlation")