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



#R5.1.1 restore stored data from own website
li="https://faculty.fordham.edu/vinod/usdata.csv"
#ourdata starts with 1969Q1 (quarter1) and ends 2006Q4
usdata=read.csv(file=li)
head(usdata,3);tail(usdata,3)#prints few lines
library(fBasics)
basicStats(usdata)
attach(usdata)

#Your registered API key is: 4d8ffce9846b8836d3f10d55a5dcd1ab
#R5.1.2  Data download from FRED.
library(fredr); fredr_set_key("4d8ffce9846b8836d3f10d55a5dcd1ab")
fredr_series_search_text("M1")
#https://fred.stlouisfed.org/docs/api/fred/series.html
devtools::install_github("keberwein/blscrapeR")
library(blscrapeR)
ids <- search_ids(keyword = c("Unemployment Rate"), periodicity_code = "q")
#Get US Money Supply M1 e dorfrom FRED
M1_dataset<-fredr(
  series_id = "M1",
  frequency = "q",
  observation_start = as.Date("1983-01-01"),
  observation_end = as.Date("2019-12-31"))
GNP_dataset<-fredr( #get GNP data
  series_id = "GNP",
  frequency = "q",
  observation_start = as.Date("1983-01-01"),
  observation_end = as.Date("2019-12-31"))
GDPDEF_dataset<-fredr(#GNP deflator inflation data
  series_id = "GDPDEF",
  frequency = "q",
  observation_start = as.Date("1983-01-01"),
  observation_end = as.Date("2019-12-31"))
#Get US unemployment data from FRED
UNRATE_dataset<-fredr(
  series_id = "UNRATE",
  frequency = "q",
  observation_start = as.Date("1983-01-01"),
  observation_end = as.Date("2019-12-31"))
#Get Import price index data from FRED
IR_dataset<-fredr(
  series_id = "IR",
  frequency = "q",
  observation_start = as.Date("1983-01-01"),
  observation_end = as.Date("2019-12-31"))



#R5.1.2 R function to convert monthly data to quarterly
#works by appending NAs and data length a multiple of 12
#If you want to average over the months set cum=F
m2qNAappend=function(x, cum=T, start.month=1,
                     end.month=12){
  y2=x
  if (start.month !=1) y2=c(rep(NA,(start.month-1)),x)
  #print(y2)
  if (end.month !=12) y2=c(y2, rep(NA,(12-end.month)))
  #print(y2)
  n=length(y2)
  #attaches NAs at the end till n is multiple of 3
  n1=n%/%3; y=y2; n2=n%%3
  if (n2>0){y=c(y2, rep(NA,n2));n1=n1+1}
  y=matrix(y,nrow=3, ncol=n1);
  #print(y)
  if (cum) quar=apply(y,2,sum, na.rm=T)
  if (cum==F) quar=apply(y,2,mean, na.rm=T)
  list(quar=quar)}
#example #x=1:18; x[3]=NA; m2qNAappend(x,start.month=2)


layout(matrix(1:4, nrow = 2, ncol = 2))
pgnp=ts(pgnp,start=c(1969,3), frequency=4)
plot(pgnp, main="US GNP Implicit Price Deflator")
ImpPr=ts(ImpPr,start=c(1969,3), frequency=4)
plot(ImpPr, main="Import Prices")
unem=ts(unem,start=c(1969,3), frequency=4)
plot(unem, main="US Unemployment rate")
realwage=ts(realwage,start=c(1969,3), frequency=4)
plot(realwage, main="US Real Wages")



#R5.1.8  Using AIC,HQ,SC, FPE criteria of VARselect command
#Needs usdata object from #R5.1.1b,2,4,5,7 in R's memory
library(vars)
VARselect(usdata, lag.max = 5, type = "const")
#$selection  #AIC(n) HQ(n) SC(n) FPE(n)
#AIC(n)  HQ(n)  SC(n) FPE(n)
#    2      2      1      2
# We choose lag=p=2 in the following VAR command
var.2c <- VAR(usdata, p = 2, type = "const")
summary(var.2c)# notation 2c for p=2 and const
plot(var.2c) #creates many plots and waits for reader to go to next
roots(var.2c)# Since root 1.0098>1 our model is unstable.


#R5.1.9  Granger-causality testing of VAR model Canadian data
library(vars); data(Canada)
#first fit the VAR model to the data choosing lag order=2
#const will include the deterministic regressor column of ones
var.2c <- VAR(Canada, p = 2, type = "const")
#output is an R object named (.2c) for 2 lags and c for const
causality(var.2c, cause = "e")


#R5.1.10 Out-of-sample forecasting of VAR models and fan-chart
#Needs usdata object from #R5.1.1b,2,4,5,7,8 in R's memory
var.2c <- VAR(usdata, p = 2, type = "const")
var.f10 <- predict(var.2c, n.ahead = 10, ci = 0.95)
names(var.f10)
layout(matrix(1:3, nrow = 3, ncol = 1))
plot(var.f10)
#fanchart(var.f10)  This is an alternative command



#R5.1.11  Impulse response function for a fitted R and plots
#Needs usdata object from #R5.1.1b,2,4,5,7,8,10 in R's memory
#irf2c=irf(var.2c) #computes all possible IRFs
colnames(usdata)  # should see following line
cnam =c("money", "realgnp",  "unem", "realwage", "pgnp", "ImpPr")
irf2c=irf(var.2c, impulse = "money",
          response = c("realgnp", "unem", "realwage", "pgnp", "ImpPr"),
          boot =F)
irf2c$irf$mo  #prints all responses to money impulse
#study the 11 by 5 matrix of responses
#length(irf2c$irf$mo)# Note two dollar symbols mo=money
nstart=seq(1,55,11)#n.ahead=10 plots 11 periods out
nend=11*1:5
# layout for 3 plots in one column
layout(matrix(1:3, nrow = 3, ncol = 1))
for ( i in 1:3){
  plot.ts(irf2c$irf$mo[nstart[i]:nend[i]],
          main="Response to Impulse in Money Supply",
          ylab=cnam[i+1])}
# layout for 2 plots in one column
layout(matrix(1:2, nrow = 2, ncol = 1))
for ( i in 4:5){
  plot.ts(irf2c$irf$mo[nstart[i]:nend[i]],
          main="Response to Impulse in Money Supply",
          ylab=cnam[i+1])}



#R5.2.1 Canonical Correlations estimation of VAR on 5
#variables of Modigliani model always available in R
rm(list=ls()) #clean out memory of R new problem
attach(LifeCycleSavings) #make names accessible
reg1=lm(sr~dpi+ddpi+pop15+pop75)
#(Intercept)         dpi        ddpi       pop15       pop75
# 28.5660865  -0.0003369   0.4096949  -0.4611931  -1.6914977
# t's   3.9      -0.36          2.1      -3.2     -1.6
cc=cancor(cbind(pop15,pop75), cbind(sr,dpi,ddpi))
# x first then y is counterintuitive to regressing y on x,
#  but this is the way it is in R
cc
alph=cc$xcoef[,1]
bet=cc$ycoef[,1]
if (bet[1]<0) {bet=-bet; alph=-alph}
rho=sqrt(cc$cor[1])
print("The fitted relation is:")
cout=c(bet,rho*alph)
round(cout,5) #round answers for printing
#      sr      dpi     ddpi    pop15    pop75
# 0.00847  0.00013  0.00417 -0.00827  0.04418



#R5.2.2  bootstrap confidence intervals (Lifecycle)
attach(LifeCycleSavings) #data into current R memory
n999=999; nx=2; ny=3; nall=nx+ny
cccoef=matrix(rep(NA, nall*n999), nrow=n999, ncol=nall) #place
#to store answers
set.seed(234)
for (i in 1:n999){
  sr2=sample(sr,replace=T)
  dpi2=sample(dpi,replace=T)
  pop152=sample(pop15,replace=T)
  ddpi2=sample(ddpi,replace=T)
  pop752=sample(pop75,replace=T)
  cc=cancor(cbind(pop152,pop752),  cbind(sr2,dpi2,ddpi2))
  # x first then y
  alph=cc$xcoef[,1]
  bet=cc$ycoef[,1]
  rho=sqrt(cc$cor[1])
  #print("The fitted relation is:")
  #if sign of weight on sr is negative, change all signs
  if (bet[1]<0) {bet=-bet; alph=-alph}
  #print(c(bet,rho*alph))
  cccoef[i,]=c(bet,rho*alph)}
for (j in 1:nall){x=sort(cccoef[,j])
cccoef[,j]=x} #sort all columns individually
print("Lower and upper Limit of a 95% confidence interval")
lolim=round(cccoef[round((n999)*0.025),],4)
uplim=round(cccoef[round((n999)*0.975),],4)
nam=c(colnames(cbind(sr2,dpi2,ddpi2)),colnames(cbind(pop152,
                                                     pop752)))
print(cbind(nam,lolim,uplim),q=F)



#R5.2.3  Joint production of Wool and Mutton
rm(list=ls()) #clean out memory of R new problem
labor=c(17703, 18552, 18814, 13490, 13390, 12710, 13100,
        13350, 14280, 14710, 13290, 12560)
cap=c(193649, 227548, 206165, 261888, 260120, 158990,
      262540, 286690, 307320, 309180, 302420, 304810)
wool=c(228091, 233309, 232258, 235807, 241284, 242177,
       239101, 243713, 259939, 166563, 261249, 248538)
mutton=c(11416, 14304, 16321, 16255, 16553, 16328, 15292,
         14495, 15528, 16239, 17536, 17171)
labor=ts(labor, frequency=1, start=c(1951,1))
cap=ts(cap, frequency=1, start=c(1951,1))
wool=ts(wool, frequency=1, start=c(1951,1))
mutton=ts(mutton, frequency=1, start=c(1951,1))
nn=10000  #used to re-scale all large data numbers
cap=cap/nn; labor=labor/nn; wool=wool/nn; mutton=mutton/nn
# Now prepare to plot data
layout(matrix(1:2, nrow = 2, ncol = 1))
plot.ts(wool, main="US Production of Wool")
plot.ts(mutton, main="US Production of Mutton")
#wait and see the first two plots
plot.ts(cap, main="Capital input in joint production")
plot.ts(labor, main="Labor input in joint production")
cc=cancor(cbind(cap,labor), cbind(wool,mutton))
# x first then y is counterintuitive to regressing y on x,
#  but this is the way it is in R
cc  #print canonical correlation results.
alph=cc$xcoef[,1]
bet=cc$ycoef[,1]
#following suggests how to fix sign problem
if (bet[1]<0) {bet=-bet; alph=-alph}
rho=sqrt(cc$cor[1])
print("The fitted relation is:")
orig=c(bet,rho*alph)
orig #original coefficients before any bootstrap
#check that we got best fitting linear combination
dewool=wool-mean(wool) #de =de-mean operation
demutt=mutton-mean(mutton)
decap=cap-mean(cap)
delab=labor-mean(labor)
cbind((orig[1]*dewool+orig[2]*demutt),
      (orig[3]*decap+orig[4]*delab))
mean((orig[1]*dewool+orig[2]*demutt)-
       (orig[3]*decap+orig[4]*delab))
#average error is zero
n999=999; nx=2; ny=2; nall=nx+ny
cccoef=matrix(rep(NA, nall*n999), nrow=n999, ncol=nall)
#place to store answers
#  Now maximum entropy bootstrap conf intervals
library(meboot)
set.seed(234) #first create 999 time series for each
wool2=meboot(x=wool, reps=n999, trim=0.10, elaps=F)
mutton2=meboot(x=mutton, reps=n999, trim=0.10, elaps=F)
cap2=meboot(x=cap, reps=n999, trim=0.10, elaps=F)
labor2=meboot(x=labor, reps=n999, trim=0.10, elaps=F)
for (i in 1:n999){ #begin loop  meboot resamples
  capz=cap2$ens[,i]; laborz=labor2$ens[,i]
  woolz=wool2$ens[,i]; muttonz=mutton2$ens[,i]
  cc=cancor(cbind(capz,laborz), cbind(woolz,muttonz))
  alph=cc$xcoef[,1]
  bet=cc$ycoef[,1]
  rho=sqrt(cc$cor[1])
  #print("The fitted relation is:")
  #if sign of weight on sr is negative, change all signs
  if (bet[1]<0) {bet=-bet; alph=-alph}
  #print(c(bet,rho*alph))
  cccoef[i,]=c(bet,rho*alph)#store the coefficients
} #end the loop using meboot resamples
#Now find confidence intervals from resamples
for (j in 1:nall){x=sort(cccoef[,j])
cccoef[,j]=x} #sort all columns individually and end j loop
print("Lower and upper Limit of a 95% confidence interval")
lolim=round(cccoef[round((n999)*0.025),],4)#lower lim at 2.5%
uplim=round(cccoef[round((n999)*0.975),],4)#upper at 97.5%
nam=c(colnames( cbind(wool,mutton)),colnames(cbind(cap,labor)))
print("Name, Lower, Original, Upper")
orig=round(orig,4)#round original coefficients for printing
print(cbind(nam,lolim,orig,uplim),q=F)#output follows
#     nam    lolim   orig    uplim
#[1,] wool   0.0127  0.0312  0.3352
#[2,] mutton -2.2689 1.7236  2.8714
#[3,] cap    -0.041  0.0119  0.0375
#[4,] labor  -1.9498 -0.8472 0.0613


#R5.2.4 Joint production of Wool and Mutton Cobb\(-\)Douglas

#specification with logs this time. Re-enter the data
labor=c(17703, 18552, 18814, 13490, 13390, 12710, 13100,
        13350, 14280, 14710, 13290, 12560)
cap=c(193649, 227548, 206165, 261888, 260120, 158990,
      262540, 286690, 307320, 309180, 302420, 304810)
wool=c(228091, 233309, 232258, 235807, 241284, 242177,
       239101, 243713, 259939, 166563, 261249, 248538)
mutton=c(11416, 14304, 16321, 16255, 16553, 16328, 15292,
         14495, 15528, 16239, 17536, 17171)

labor=log(labor)#logs needed for the Cobb\(-\)Douglas model

cap=log(cap)
wool=log(wool)
mutton=log(mutton)
#  define data as time series
labor=ts(labor, frequency=1, start=c(1951,1))
cap=ts(cap, frequency=1, start=c(1951,1))
wool=ts(wool, frequency=1, start=c(1951,1))
mutton=ts(mutton, frequency=1, start=c(1951,1))
tim=51:62  #index of technological change
cc=cancor(cbind(cap,labor,(cap*labor)),cbind(wool,mutton))
# x first then y is counterintuitive to regressing y on x,
#  but this is the way it is in R
cc
alph=cc$xcoef[,1]
bet=cc$ycoef[,1]
if (bet[1]<0) {bet=-bet; alph=-alph}
rho=sqrt(cc$cor[1])
print("The fitted relation is:")
orig=c(bet,rho*alph)
orig
# copying the code on meboot from #R5.2.3 gives conf intervals as:
#       nam    lolim   orig    uplim
#wool   wool   0.2543  0.6624  8.1705
#mutton mutton -3.2939 2.4864  4.6995
#cap    cap    -0.9745 0.2625  0.8578
#labor  labor  -2.7674 -1.3677 0.2599

# Thus, Cobb\(-\)Douglas intervals similar to those without logs

#extra regressor for technology makes coeff for cap and labor
#both negative.  Details omitted for brevity.



#R5.2.5: CCA for generalized correlations
#R5.2.2 Life cycle data should be in memory
mtx=cbind(pop15,pop75,sr,dpi,ddpi)
options(np.messages=FALSE);Rstar=gmcmtx0(mtx)
canonRho(Rstar,nLHS=2, ridg=c(.1,.1),verbo=TRUE)


#R5.3.1 Sophisticated cointegration tests of urca package
#ablrtest =Likelihood ratio test, restrictions on alpha, beta
#Ukpppuip = Dataset for the United Kingdom: ppp and uip
#alphaols =OLS, VECM=vector error correction weighting matrix
#alrtest =Likelihood ratio test for restrictions on alpha
#bh5lrtest =Likelihood ratio test, restrict partly known beta
#blrtest =Likelihood ratio test for restrictions on beta
#ca.jo =Johansen Procedure for VAR, reports trace
#ca.po =Phillips \& Ouliaris Cointegration Test
#cajolst =Test Cointegrating Rank, Level Shift, Unknown time
#cajools =OLS regression of VECM
#lttest =Likelihood ratio test for no linear trend in VAR
#plotres =Graphical inspection of VECM residuals
#r.df =Augmented-Dickey-Fuller Unit Root Test



#R5.4.1 SVAR example Mpyr data
library(Ecdat);data(Mpyr);attach(data.frame(Mpyr))
library(vars)
v1=VAR(cbind(p,y),exogen=cbind(m,r),p=2,ic="AIC")
roots(v1)  #stable because all roots <1
causality(v1, cause = "p") # p does not cause y
BQ(v1)