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


#R8.3.1.1  Data Download Monthly from March 1959
lin="https://faculty.fordham.edu/vinod/ch8data.txt"
ch8data=read.table(file=lin,header=T)
head(ch8data,3)
attach(ch8data) #data names now accessible
rcnd=ts(rcnd,start=c(1959,3), frequency=12)
popn=ts(popn, start=c(1959,3), frequency=12)
pccnd=ts(pccnd, start=c(1959,3), frequency=12)


#R8.3.1.2 Creating a pdf file of data plots
#pdf(file="c:/temp/plotCY.pdf", width=6.25, height=6,
#family="Times", pointsize=12)
layout(matrix(1:3, nrow = 3, ncol = 1))
plot(rcnd, main="Real personal consumption expenditures
non-durable goods (billions, chained index, base=2000)",
     xlab="Time (monthly)")
ma="Total Population, thousands, All Ages + Armed Forces
Overseas"
plot(popn, main=ma,xlab="Time  (monthly)")
ma="per capita consumption of non-durables"
plot(pccnd, main=ma, xlab="Time (monthly)")
#dev.off()  #this command needed to stop the pdf file
ar(pccnd) #fit AR(1)


#R8.3.2 OLS regression and tests for autocorrelation and
#heteroscedasticity
library(dyn)
reg1=dyn$lm(pccnd~lag(pccnd,-1)) #$
#ar(consu, method="ols")#not clear different results
summary(reg1)
#Breusch-Godfrey test for serial correlation of errors
library(lmtest)
for (ord in 1:4){
  print(bgtest(reg1, ord))}
#library(xtable) to export the table to LATEX
#xtable(summary(reg1))
library(car)
bptest(reg1)# studentized Breusch-Pagan test



#R8.3.3.1 Function to compute concentrated likelihood
fn<-function(xx, pccnd) {#name of function is fn
  psi=xx[1]
  alpha=xx[2]^2  #Square the second coefficient of answer
  #vector to get a positive alpha>0
  n=length(pccnd); constant=1; psi2=psi^(-1/alpha)
  alp2=-2*alpha-2; alp3=2*alpha+2
  tby2=(n-1)/2 #t1 term 1 ct - psi-1/alph
  s1=rep(NA,n-1); s2=rep(NA,n-1)
  for(i in 2:n)
    s1[i-1]={((pccnd[i]-psi2*pccnd[i-1])^2)*(pccnd[i-1]
                                             ^alp2)}
  for(m in 2:n)
  {s2[m-1]=log(pccnd[m-1])^alp3}
  s3<-constant-tby2*log(sum(s1))-(1/2)*sum(s2)
  return(s3) }#function fn ends here



#R8.3.3.2 Actual Optimization Calls and Printing of
#Results
oNM=optim(c(0.99,0.1),fn, method="Nelder-Mead",
          pccnd=pccnd)
oNM; alp=oNM$par[2]^2
psi=oNM$par[1]
print(c("alpha=",alp, "gamma=", psi),q=F)
alp=round(alp,8); psi=round(psi,8)
print(c("Nelder-Mead alpha=",alp, "gamma=", psi),q=F)
oBF=optim(c(0.99,0.1),fn, method="BFGS",pccnd=pccnd)
alp=oBF$par[2]^2; psi=oBF$par[1]
print(c("alpha=",alp, "gamma=", psi),q=F)
alp=round(alp,8); psi=round(psi,8)
print(c("BFGS alpha=",alp, "gamma=", psi),q=F)
#USE BHHH optimization algorithm
library(micEcon);library(maxLik)
mB=maxBHHH(fn,start=c(0.99,0.1),pccnd=pccnd,
           iterlim=300, grad=NULL)  #maxBHHH failed here
#mB; alp=mB$estimate[2]^2
#psi=mB$estimate[1]
#print(c("alpha=",alp, "gamma=", psi),q=F)
#alp=round(alp,8); psi=round(psi,8)
#print(c("BHHH alpha=",alp, "gamma=", psi),q=F)


#R8.3.3.3 Quarterly data estimation of parameters
pccnd2=m2q(pccnd)$quar  #convert to quarterly series
# you need the "$quar" in the previous command to get the
#series out
oNM2=optim(c(0.99,0.1),fn, method="Nelder-Mead",
           pccnd=pccnd2)
oNM2; alp=oNM2$par[2]^2; psi=oNM2$par[1]
print(c("oNM2 alpha=",alp, "gamma=", psi),q=F)
alp=round(alp,8); psi=round(psi,8)
print(c("oNM2 alpha=",alp, "gamma=", psi),q=F)
oBF2=optim(c(0.99,0.1),fn, method="BFGS",pccnd=pccnd2)
alp=oBF2$par[2]^2;psi=oBF2$par[1]
print(c("oBF2 alpha=",alp, "gamma=", psi),q=F)
alp=round(alp,8); psi=round(psi,8)
print(c("oBF2 alpha=",alp, "gamma=", psi),q=F)
#USE BHHH optimization algorithm
library(micEcon)
mB2=maxBHHH(fn,theta=c(0.99,0.1),pccnd=pccnd2,iterlim=150)
mB2; alp=mB2$estimate[2]^2; psi=mB2$estimate[1]
print(c("BHHH2 alpha=",alp, "gamma=", psi),q=F)
alp=round(alp,8); psi=round(psi,8)
print(c("BHHH2 alpha=",alp, "gamma=", psi),q=F)



#R8.4.1 non-parametric kernel regression
library(np); options(np.messages=FALSE)
bw <- npregbw(formula=pccnd[2:n]~pccnd[1:(n-1)], tol=.1,
              ftol=.1)
# Kernel Regression Bandwidth Selection with Mixed data
n=length(pccnd)
# Now take these bandwidths and fit the model and
# gradients
model <- npreg(bws = bw, gradients = TRUE,
               residuals=T)
# Kernel Regression with Mixed Datatypes
myresi=model$resid
mstar=myresi-mean(myresi) #unanticipated pccndmption
head(mstar); tail(mstar)
head(model$eval)  #evaluation points
tail(model$eval)  #evaluation points
head(model$grad)
tail(model$grad)#this evaluates gradients
matplot(cbind(model$eval,model$mean), typ="l",
        main="C(t-1)and fitted values",
        ylab="data and fitted values")
matplot(cbind(model$merr,model$grad), typ="l",
        main="standard errors and gradients",
        ylab=" SE and gradients")
model$R2  # R square coeff of determination
head(cbind(pccnd[2:n], (pccnd[2:n]-model$resid))) #$
plot.ts(mstar, main="Unanticipated consumption")
start(pccnd);end(pccnd) #to get the time series
#properties
unanti=ts(mstar,start=c(1959,4),frequency=12)
end(unanti) #how far does the data go?
cor(rdispi[2:n],unanti) #simple corr requires care in
#making
#data comparable, ccf is smart, makes vectors comparable
#[1] 0.02533918 low correlations with real disposable
#income. ccf=cross correlation function.
ccf(unanti,prime, lag.max=2,plot=F)#lags of first series
# -0.052  -0.057  -0.055  -0.059  -0.061
ccf(unanti,rdispi, lag.max=2,plot=F) #lags of first
#series
#0.028   0.027   0.025   0.018   0.020


q=seq(from=1,to=10,by=.1)
p=50-3*q
plot(q,p,typ="l",main="Demand Curve p=50-3q")
p1=29; p2=35; q1=7; q2=5
lines(x=c(0,q1),y=c(p1,p1))
lines(x=c(0,q2),y=c(p2,p2))
pdiff=abs(p2-p1)
qmin=min(q1,q2)
qmax=max(q1,q2)
lines(x=c(qmin,qmin),y=c(p1,p2))
tribase=qmax-qmin#base of triangle
CS=(qmin*pdiff)+0.5*(pdiff*tribase);CS



#R8.A.1 download data as in \#R8.3.1.1.
head(cbind(percapc,percapi))
#get ready of Granger causality study
library(vars)
myy=data.frame(cbind(percapc,percapi))
var.1c <- VAR(myy, p = 1, type = "const")
#output is an R object named (var.1c) where
#we have 1 lag and c is for const (intercept)
causality(var.1c, cause = "percapi")
causality(var.1c, cause = "percapc")
# results show that consumption Granger-causes income,
#not
#vice versa
myy=data.frame(cbind(pccnd,percapi))
var.1c <- VAR(myy, p = 1, type = "const")
#output is an R object named (.1c) for 1 lag and c for
#const
causality(var.1c, cause = "percapi")
causality(var.1c, cause = "pccnd")
ccf(unanti,prime, lag.max=2,plot=F)
# -0.052  -0.057  -0.055  -0.059  -0.061
ccf(unanti,rdispi, lag.max=2,plot=F) #lags of first
#series
#0.028   0.027   0.025   0.018   0.020


#R8.A.2 function for converting monthly to quarterly
#series
m2q=function(x, cum=T){
  # Input= x monthly series, OUTPUT= quar as quarterly
  #version
  #of x
  y2=x; n=length(x)
  start.year=start(x)[1]
  new.year=start.year
  start.month=start(x)[2]
  end.month=end(x)[2]
  end.year=end(x)[1]
  if (start.month==1) new.start=1
  if (start.month==2) {new.start=2
  #delete Feb and March first two data points
  y2=window(x, start=c(new.year,4))}
  if (start.month==3) {new.start=2
  #delete March first  1 data point
  y2=window(x, start=c(new.year,4))}
  if (start.month==4) new.start=2
  if (start.month==5) {new.start=3
  #delete May and June 2 data point
  y2=window(x, start=c(new.year,7))}
  if (start.month==6) {new.start=3
  #delete June  1 data point
  y2=window(x, start=c(new.year,7))}
  if (start.month==7) new.start=3
  if (start.month==8) {new.start=4
  #delete Aug and Sept 2 data points
  y2=window(x, start=c(new.year,10))}
  if (start.month==9) {new.start=4
  #delete Sept 1 data point
  y2=window(x, start=c(new.year,10))}
  if (start.month==10) new.start=4
  if (start.month==11) {new.start=1
  #delete NOV and DEC first  2 data points
  new.year=start.year+1
  y2=window(x, start=c(new.year,1))}
  if (start.month==12) {new.start=1
  #delete DEC first  1 data point
  new.year=start.year+1
  y2=window(x, start=c(new.year,1))}
  eyy2=end(y2)[1]
  emy2=end(y2)[2]
  nn=length(y2)
  if((nn%%3)!=0) y2=window(y2, end=c(eyy2-1,12))
  print(c("starting value for selected monthly series",
          start(y2)),q=F)
  print(c("ending value for selected monthly series",
          end(y2)),q=F)
  n=length(y2)
  y=matrix(y2,nrow=3);
  #print(y) #for details remove the # at the start of
  #this line
  if (cum) quar=apply(y,2,sum, na.rm=F)
  if (cum==F) quar=apply(y,2,mean, na.rm=F)
  quar=ts(quar, start=c(new.year,new.start), frequency=4)
  print(c("starting value for quarterly series",
          start(quar)),q=F)
  print(c("ending value for quarterly series",
          end(quar)),q=F)
  list(quar=quar)} #function m2q ends here
#example
#x=ts(1:36, start=c(1959,2), frequency=12); x[2]=NA; x;
#m2q(x)


#R8.A.3 Fit VAR model to quarterly version of the data
lin="https://faculty.fordham.edu/vinod/ch8quar.txt"
ch8quar=read.table(file=lin,header=T)
head(ch8quar,3)
attach(ch8quar) #pc=per-capita q=quarterly
myy=data.frame(cbind(pcCq,pcYq))
library(vars)
var.1cq <- VAR(myy, p = 1, type = "const")
#output is an R object named (.1c) for 1 lag and c for
#const
causality(var.1cq, cause = "pcCq")
causality(var.1cq, cause = "pcYq")



#R8.A.4 Target consumption construction
lin="https://faculty.fordham.edu/vinod/ch8data.txt"
ch8data=read.table(file=lin,header=T)
head(ch8data,3); attach(ch8data) #
library(dyn);n=length(infl);n; tim=1:n
reg2=lm(percapc[2:n]~percapc[1:(n-1)]+percapi[2:n]+
          infl[2:n]+tim[2:n])
summary(reg2)
head(reg2$coef[1]+as.matrix(reg2$model[,2:5]) %*%
       as.vector
     (reg2$coef[2:5]))
head(fitted(reg2))
#above two agree so we can proceed to ignore the
#time regressor with a negative coeff
head(reg2$coef[1]+as.matrix(reg2$model[,2:4]) %*%
       as.vector
     (reg2$coef[2:4]))
#Now define C*
cstar=(reg2$coef[1]+as.matrix(reg2$model[,2:4]) %*%
         as.vector
       (reg2$coef[2:4]))
#note cstar > consumption as required by theory
#Now fit the decision equation of Wiener-Hopf theory
reg3=lm(percapi[2:n]~percapc[1: (n-1)]+cstar);
#summary(reg3)
q1=reg3$coef[2]
q2=reg3$coef[3]
#NOw fit the habit equation
reghabit=lm(percapc[2:n]~percapi[2:n]+percapc[1:(n-1)])
b1=reghabit$coef[1]
b2=reghabit$coef[2]
b3=reghabit$coef[3]
#compute variance of income.
vinc=var(percapi[2:n])
#compute variance of C - C*
varccstar=var(percapc[2:n]-cstar)
mu1=b2*(q1*b2+b3)/ (q2*b3)
mu2=-2* (mean(cstar)-mean(percapc[1: (n-1)]) )/q1
mu1;mu2 #$