This is the program file for the article "Rising Natural Gas Prices and Real Economc Activity" 
by Kevin L. Kliesen in the November/December 2006 issue of the Review.


####################################### Functions ##################################################

# Save these functions in a separate text file named R_functions.txt:

### Function to calculate lags for a vector ###
lags<- function(obj,lag=1) {
   v = c(rep(NA,lag),obj)
   return( v[1:length(obj)] )
}
#---------------# End of function #---------------#

### Function to calculate changes for a matrix ###
delta <- function(x, n=1, type=1) {

x = as.matrix(x)
nrows = dim(x)[1]; n0=nrows; n1=n0-n;
ncols = dim(x)[2]
tmp = matrix(NA, nrow=nrows, ncol=ncols)

for(i in 1:ncols) {
   if(type==1) tmp[,i] = c( rep(NA,n),     x[(1+n):n0,i]/x[1:n1,i]-1 )  # discrete price changes
   if(type==2) tmp[,i] = c( rep(NA,n), log(x[(1+n):n0,i]/x[1:n1,i])  )  # continuous volume changes
   if(type!=1 & type!=2)
      stop("Invalid compounding type \n", "type = 1 for discrete, 2 for continuous")
}
return(tmp)
}
#---------------# End of function #---------------#

########################################## Table 4 ##################################################

# Program written:
# in: R1.9.1
#
# This program runs regressions similar to Hamilton 2003's regression 3.2, except the
# energy data are just lags of the first log difference (not Hamilton transformed) over
# all the 3 digit IP industries and total IP.  It runs the regressions 12 times for lags
# 1; 1:2, 1:3, ..., 1:12.  Then it chooses the optimal lag structure as the one with the
# lowest AIC statistic.  Finanlly, it performes an F-test on the regression with the
# lowest AIC to test if all the energy lags equal zero.

# Copy the yellow highlighted data from "Data.xls" (sheet: Table 4)
# and paste it in a text file named "3Digit.Reg.Data.txt" 
# The ouput file is "AllX.Reg.Output.txt" where X = Oil or Gas.

### Choose which PPI you want to run regressions with ###
which.PPI = 1  ### 1 = gas, 2 = oil

# Change directory as necessary:
setwd("H:/Kliesen/Published Papers/Review/2006/NG Paper/NG Stuff/Test")

### Function to calculate lags for a vector ###
lags<- function(obj,lag=1) { v <- c(rep(NA,lag),obj); v[1:length(obj)] }

### Function to calculate changes for a matrix ###
delta <- function(x, n=1, type=1) {

x = as.matrix(x)
nrows = dim(x)[1]; n0=nrows; n1=n0-n;
ncols = dim(x)[2]
tmp = matrix(NA, nrow=nrows, ncol=ncols)

for(i in 1:ncols) {
   if(type==1) tmp[,i] = c( rep(NA,n),     x[(1+n):n0,i]/x[1:n1,i]-1 )  # discrete price changes
   if(type==2) tmp[,i] = c( rep(NA,n), log(x[(1+n):n0,i]/x[1:n1,i])  )  # continuous volume changes
   if(type!=1 & type!=2)
      stop("Invalid compounding type \n", "type = 1 for discrete, 2 for continuous")
}
return(tmp)
}
################ End of functions #################

### Read in data, and remove any "#N/A" rows ###
# Columns of data #
# 1: Date
# 2: PPI: Gas Fuels (NSA, 1982=100)
# 3: PPI: Crude Petroleum [Domestic Production] (NSA, 1982=100)
# 4: IP: Manufacturing (SIC) (SA, 1997=100) 
# 5: IP: Food (SA, 1997=100) 
# 6: IP: Beverages and Tobacco Products (SA, 1997=100) 
# 7: IP: Textile Mills (SA, 1997=100) 
# 8: IP: Textile Product Mills (SA, 1997=100) 
# 9: IP: Apparel (SA, 1997=100) 
#10: IP: Leather and Allied Products (SA, 1997=100) 
#11: IP: Wood Products (SA, 1997=100) 
#12: IP: Paper and Paper Products (SA, 1997=100) 
#13: IP: Printing and Related Support Activities (SA, 1997=100) 
#14: IP: Petroleum and Coal Products (SA, 1997=100) 
#15: IP: Chemicals (SA, 1997=100) 
#16: IP: Plastics and Rubber Products (SA, 1997=100) 
#17: IP: Nonmetallic Mineral Products   (SA, 1997=100) 
#18: IP: Primary Metals (SA, 1997=100) 
#19: IP: Fabricated Metal Products (SA, 1997=100) 
#20: IP: Machinery (SA, 1997=100) 
#21: IP: Computer and Electronic Products (SA, 1997=100) 
#22: IP: Electrical Equipment, Appliances, and Components (SA, 1997=100) 
#23: IP: Transportation Equipment (SA, 1997=100) 
#24: IP: Furniture and Related Products (SA, 1997=100) 
#25: IP: Miscellaneous Manufacturing (SA, 1997=100) 

xcols = 24
x = matrix( scan( "3Digit.Reg.Data.txt", skip=3, what="" ), byrow=T, ncol=xcols+1 )
x = matrix( as.numeric( na.omit( ifelse(x=="#N/A", NA, x) ) ), ncol=xcols+1 )

### Remove non-chosen PPI ###
if(which.PPI==1) x = x[,-3]
if(which.PPI==2) x = x[,-2]

# Take log first differences of columns 2:xcols
x[,2:xcols] = delta( x[,2:xcols], type=2, n=1 )	
colnames(x) = c("Date", "dln.PPI", "dln.IP", "dln.Food", "dln.Bev.Tobacco", "dln.Textile.Mills",
	"dln.Textile.Prod", "dln.Apparel", "dln.Leather", "dln.Wood", "dln.Paper", "dln.Printing",
	"dln.Petroleum.Coal", "dln.Chem", "dln.Plastic.Rubber", "dln.Nonmetal.Mineral",
	"dln.Primary.Metal", "dln.Fab.Metal", "dln.Machinery", "dln.Computer", "dln.Elec.Equip",
	"dln.Transport", "dln.Furniture", "dln.Misc.Manufact")
x = data.frame(x)

### Run regression(s) #5 ###
# Create 12 lags of PPI
PPIlags = NULL
for(i in 1:12)   PPIlags = cbind( PPIlags, lags(x[,2], i) )
colnames(PPIlags) = paste("dln.PPI.",1:12, sep="")

# Create 12 lags of each measure of IP
IPcolnams = NULL
IPlags = NULL

for(j in 3:xcols) {
   for(i in 1:12) {
      IPlags = cbind( IPlags, lags(x[,j], i) )
      IPcolnams = c( IPcolnams, paste(colnames(x)[j],".",i,sep="") )
   }
}
colnames(IPlags) = IPcolnams

xlags = na.omit( cbind( x, IPlags, PPIlags ) )
colnames(xlags) = c( colnames(x), colnames(IPlags), colnames(PPIlags) )

### Begin Loop over industries ###
for(i in 3:xcols) {

### Set up some Vectors ###
seq = (25+12*(i-3)):(36+12*(i-3))

ip.lags = names(xlags[seq[1:12]])
ppilags = names(xlags[(ncol(xlags)-11):ncol(xlags)])
ip.vars = NULL
ppivars = NULL
for(k in 1:12) {
ip.vars = cbind( ip.vars, gsub(","," +",toString(ip.lags[1:k])) )
ppivars = cbind( ppivars, gsub(","," +",toString(ppilags[1:k])) )
}

# Create Formulas
for(h in 1:12) {
assign(paste("f",h,sep=""),     as.formula( paste(names(xlags[i]),"~",ip.vars[h],"+", ppivars[h]) ))
assign(paste("f",h,"R",sep=""), as.formula( paste(names(xlags[i]),"~",ip.vars[h]) ))
}
# End Create Formulas

# Run Regressions
r51  = lm( f1 , data=xlags )
r52  = lm( f2 , data=xlags )
r53  = lm( f3 , data=xlags )
r54  = lm( f4 , data=xlags )
r55  = lm( f5 , data=xlags )
r56  = lm( f6 , data=xlags )
r57  = lm( f7 , data=xlags )
r58  = lm( f8 , data=xlags )
r59  = lm( f9 , data=xlags )
r510 = lm( f10, data=xlags )
r511 = lm( f11, data=xlags )
r512 = lm( f12, data=xlags )

### Find lowest AIC statistic ###
aic = AIC( r51, r52, r53, r54, r55, r56, r57, r58, r59, r510, r511, r512 )
min.aic = aic[match(min(aic[,2]),aic[,2]),]

### Get Unrestricted and Restricted Regressions for lowest AIC Statistic ###

Uformula = get( sub("r5","f",rownames(min.aic)) )
Uregress = get( rownames(min.aic) )
Uformula = gsub( ","," +",sub(","," =",toString(all.vars(Uformula))) )

Rformula = get( paste(sub("r5","f",rownames(min.aic)), "R", sep="") )
Rregress = lm( Rformula, data=xlags )
Rformula = gsub( ","," +",sub(","," =",toString(all.vars(Rformula))) )

### Do F-Test on null that all energy lags = 0 ###
Ftest = anova(Uregress, Rregress)

### Output Regression Summary Results to file ###
if(which.PPI==1) outfile = "AllGas.Reg.Output.txt"
if(which.PPI==2) outfile = "AllOil.Reg.Output.txt"
sink(outfile, append=T)  # Open connection to file

cat("\n### Unrestricted Regression:\n### ", Uformula,
    "\n### Restricted Regression:\n### ", Rformula,
    "\n### Unrestricted Output:\n")
print(summary(Uregress))
cat("\n### F-test Results:\n")
print(anova(Uregress,Rregress))
cat("\n\n")

sink()  # Close connection to file

}
### End Loop over Industries ###

############################################### Tables 5 & 6 ###############################################

# Program written:
# in: R1.9.1
#
# This program runs regressions similar to Hamilton 2003's regression 3.2, for Hamilton
# transformed natural gas and oil prices over all the 3 digit IP industries and total IP.
# It runs the regressions 12 times for lags 1; 1:2, 1:3, ..., 1:12.  Then it chooses the
# optimal lag structure as the one with the lowest AIC statistic.  Finally, it performes
# an F-test on the regression with the lowest AIC to test if all the energy lags equal zero.

# Copy the yellow highlighted data from "Data.xls" (sheet: Table 5_6)
# and paste it in a text file named "3Digit.Hamilton.Reg.Data.txt" 
# The ouput files are "Hamilton.XNG.Output.txt" & "Hamilton.XOil.Output.txt"
# where X equals 1 or 3 for the 12 month or 36 month Hamilton transformed data, respectively.

### Choose which PPI you want to run regressions with ###
which.PPI = 2  ### 1 = gas, 2 = oil

setwd("H:/Kliesen/Published Papers/Review/2006/NG Paper/NG Stuff/Test")

### Function to calculate lags for a vector ###
lags<- function(obj,lag=1) { v <- c(rep(NA,lag),obj); v[1:length(obj)] }

### Function to calculate changes for a matrix ###
delta <- function(x, n=1, type=1) {

x = as.matrix(x)
nrows = dim(x)[1]; n0=nrows; n1=n0-n;
ncols = dim(x)[2]
tmp = matrix(NA, nrow=nrows, ncol=ncols)

for(i in 1:ncols) {
   if(type==1) tmp[,i] = c( rep(NA,n),     x[(1+n):n0,i]/x[1:n1,i]-1 )  # discrete price changes
   if(type==2) tmp[,i] = c( rep(NA,n), log(x[(1+n):n0,i]/x[1:n1,i])  )  # continuous volume changes
   if(type!=1 & type!=2)
      stop("Invalid compounding type \n", "type = 1 for discrete, 2 for continuous")
}
return(tmp)
}
################ End of functions #################

### Read in data, and remove any "#N/A" rows ###
# Columns of data #
# 1: Date
# 2: 1 year (12 month) Hamilton Trans. PPI: Gas Fuels (NSA, 1982=100)
# 3: 1 year (12 month) Hamilton Trans. PPI: Crude Petroleum [Domestic Production] (NSA, 1982=100)
# 4: 3 year (36 month) Hamilton Trans. PPI: Gas Fuels (NSA, 1982=100)
# 5: 3 year (36 month) Hamilton Trans. PPI: Crude Petroleum [Domestic Production] (NSA, 1982=100)
# 6: IP: Manufacturing (SIC) (SA, 1997=100) 
# 7: IP: Food (SA, 1997=100) 
# 8: IP: Beverages and Tobacco Products (SA, 1997=100) 
# 9: IP: Textile Mills (SA, 1997=100) 
#10: IP: Textile Product Mills (SA, 1997=100) 
#11: IP: Apparel (SA, 1997=100) 
#12: IP: Leather and Allied Products (SA, 1997=100) 
#13: IP: Wood Products (SA, 1997=100) 
#14: IP: Paper and Paper Products (SA, 1997=100) 
#15: IP: Printing and Related Support Activities (SA, 1997=100) 
#16: IP: Petroleum and Coal Products (SA, 1997=100) 
#17: IP: Chemicals (SA, 1997=100) 
#18: IP: Plastics and Rubber Products (SA, 1997=100) 
#19: IP: Nonmetallic Mineral Products   (SA, 1997=100) 
#20: IP: Primary Metals (SA, 1997=100) 
#21: IP: Fabricated Metal Products (SA, 1997=100) 
#22: IP: Machinery (SA, 1997=100) 
#23: IP: Computer and Electronic Products (SA, 1997=100) 
#24: IP: Electrical Equipment, Appliances, and Components (SA, 1997=100) 
#25: IP: Transportation Equipment (SA, 1997=100) 
#26: IP: Furniture and Related Products (SA, 1997=100) 
#27: IP: Miscellaneous Manufacturing (SA, 1997=100) 

xcols = 27
x = matrix( scan( "3Digit.Hamilton.Reg.Data.txt", skip=3, what="" ), byrow=T, ncol=xcols )
x = matrix( as.numeric( na.omit( ifelse(x=="#N/A", NA, x) ) ), ncol=xcols )

### Remove non-chosen PPI ###
if(which.PPI==1) x = x[,-c(3,5)]  # Run Gas, remove Oil
if(which.PPI==2) x = x[,-c(2,4)]  # Run Oil, remove Gas
xcols = ncol(x)

# Take log first differences of columns 4:xcols
x[,4:xcols] = 100 * delta( x[,4:xcols], type=2, n=1 )	
colnames(x) = c("Date", "dln.PPI.H1", "dln.PPI.H3", "dln.IP", "dln.Food", "dln.Bev.Tobacco",
        "dln.Textile.Mills", "dln.Textile.Prod", "dln.Apparel", "dln.Leather", "dln.Wood",
	"dln.Paper", "dln.Printing", "dln.Petroleum.Coal", "dln.Chem", "dln.Plastic.Rubber",
	"dln.Nonmetal.Mineral", "dln.Primary.Metal", "dln.Fab.Metal", "dln.Machinery",
	"dln.Computer", "dln.Elec.Equip", "dln.Transport", "dln.Furniture", "dln.Misc.Manufact")
x = data.frame(x)

for(m in 1:2) {  ### Begin Loop over Hamilton Transformed Data ###

### Run regression(s) #5 ###
# Create 12 lags of PPI
PPIlags = NULL
for(i in 1:12)   PPIlags = cbind( PPIlags, lags(x[,m+1], i) )
colnames(PPIlags) = paste("dln.PPI",m^2,".",1:12, sep="")

# Create 12 lags of each measure of IP
IPcolnams = NULL
IPlags = NULL

for(j in 4:xcols) {
   for(i in 1:12) {
      IPlags = cbind( IPlags, lags(x[,j], i) )
      IPcolnams = c( IPcolnams, paste(colnames(x)[j],".",i,sep="") )
   }
}
colnames(IPlags) = IPcolnams

xlags = na.omit( cbind( x, IPlags, PPIlags ) )
colnames(xlags) = c( colnames(x), colnames(IPlags), colnames(PPIlags) )

### Begin Loop over industries ###
for(i in 4:xcols) {

### Set up some Vectors ###
seq = ( (xcols + 1) + 12 * (i-4) ):( (xcols + 13) + 12 * (i-4) )

ip.lags = names(xlags[seq[1:12]])
ppilags = names(xlags[(ncol(xlags)-11):ncol(xlags)])
ip.vars = NULL
ppivars = NULL
for(k in 1:12) {
ip.vars = cbind( ip.vars, gsub(","," +",toString(ip.lags[1:k])) )
ppivars = cbind( ppivars, gsub(","," +",toString(ppilags[1:k])) )
}

# Create Formulas
for(h in 1:12) {
assign(paste("f",h,sep=""),     as.formula( paste(names(xlags[i]),"~",ip.vars[h],"+", ppivars[h]) ))
assign(paste("f",h,"R",sep=""), as.formula( paste(names(xlags[i]),"~",ip.vars[h]) ))
}
# End Create Formulas

# Run Regressions
r51  = lm( f1 , data=xlags )
r52  = lm( f2 , data=xlags )
r53  = lm( f3 , data=xlags )
r54  = lm( f4 , data=xlags )
r55  = lm( f5 , data=xlags )
r56  = lm( f6 , data=xlags )
r57  = lm( f7 , data=xlags )
r58  = lm( f8 , data=xlags )
r59  = lm( f9 , data=xlags )
r510 = lm( f10, data=xlags )
r511 = lm( f11, data=xlags )
r512 = lm( f12, data=xlags )

### Find lowest AIC statistic ###
aic = AIC( r51, r52, r53, r54, r55, r56, r57, r58, r59, r510, r511, r512 )
min.aic = aic[match(min(aic[,2]),aic[,2]),]

### Get Unrestricted and Restricted Regressions for lowest AIC Statistic ###

Uformula = get( sub("r5","f",rownames(min.aic)) )
Uregress = get( rownames(min.aic) )
Uformula = gsub( ","," +",sub(","," =",toString(all.vars(Uformula))) )

Rformula = get( paste(sub("r5","f",rownames(min.aic)), "R", sep="") )
Rregress = lm( Rformula, data=xlags )
Rformula = gsub( ","," +",sub(","," =",toString(all.vars(Rformula))) )

### Do F-Test on null that all energy lags = 0 ###
Ftest = anova(Uregress, Rregress)

### Output Regression Summary Results to file ###
if(which.PPI==1) outfile = paste( "Hamilton.",m^2, "NG.Output.txt", sep="")
if(which.PPI==2) outfile = paste( "Hamilton.",m^2,"Oil.Output.txt", sep="")
sink(outfile, append=T)  # Open connection to file

cat("\n### Unrestricted Regression:\n### ", Uformula,
    "\n### Restricted Regression:\n### ", Rformula,
    "\n### Unrestricted Output:\n")
print(summary(Uregress))
cat("\n### F-test Results:\n")
print(anova(Uregress,Rregress))
cat("\n\n")

sink()  # Close connection to file

### Output Regression Lags, F-statistic, & P-value to table ###
ip.sector = substr( unlist(strsplit(Uformula," ="))[1], 5, 25 )
aic.lags = as.numeric( substr(row.names(min.aic),3,nchar(row.names(min.aic))) )
fstat = anova(Uregress,Rregress)[2,5]
pvalue = anova(Uregress,Rregress)[2,6]

if(which.PPI==1) name = "NG";  if(which.PPI==2) name = "Oil";
outtable = paste( "Hamilton.",name,".Table.Output.txt", sep="")
sink(outtable, append=T)
cat( m^2, ip.sector, aic.lags, fstat, pvalue, "\n" )
sink()

}
### End Loop over Industries ###

}
### End Loop over Hamilton Transformed Data ###

################################################### Table 7 ############################################

# Program written:
# in: R1.9.1
#
# This program replicates the regressions 3.2 & 3.8 from Hamilton 2003, and then runs them
# over Gas PPI data.  It also calculates an F-Test that tests if the coefficients on the
# lagged gas prices are all equal to zero for both regressions.

# Copy the yellow highlighted data from "Data.xls" (sheet: Table 7) 
# and save it in a text file named "Hamilton JOE 2003 Data.txt."
# The output file is "Output_Compare Hamilton Oil V NG.txt."

# Change these directories as necessary:
setwd("H:/Kliesen/Published Papers/Review/2006/NG Paper/NG Stuff/Test")
source("H:/Kliesen/Published Papers/Review/2006/NG Paper/NG Stuff/Test/R_functions.txt")
outfile = "Output_Compare Hamilton Oil V NG.txt"

### Read in data, and replace "#N/A" with NA ###
xcols = 11
x = matrix(scan("Hamilton JOE 2003 Data.txt", what="", sep="\t"), byrow=T, ncol=xcols )
#x = matrix( as.numeric( na.omit( ifelse(x=="#N/A", NA, x) ) ), ncol=xcols )
x = matrix( as.numeric( ifelse(x=="#N/A", NA, x) ), ncol=xcols )
colnames(x) = c("Date","rGDP","PPI_Oil","PPI_Gas","dln.rGDP","dln.PPI_Oil","dln.PPI_Gas",
                "H1.PPI_Oil","H1.PPI_Gas","H3.PPI_Oil","H3.PPI_Gas")

x[,"dln.rGDP"] = 100 * delta(x[,"rGDP"], type=2)
x[,"dln.PPI_Oil"] = 100 * delta(x[,"PPI_Oil"], type=2)
x[,"dln.PPI_Gas"] = 100 * delta(x[,"PPI_Gas"], type=2)

oil.data = cbind(x[,"Date"], x[,"dln.rGDP"], x[,"dln.PPI_Oil"], x[,"H1.PPI_Oil"], x[,"H3.PPI_Oil"])
colnames(oil.data) = colnames(x)[c(1,5,6,8,10)]
gas.data = cbind(x[,"Date"], x[,"dln.rGDP"], x[,"dln.PPI_Gas"], x[,"H1.PPI_Gas"], x[,"H3.PPI_Gas"])
colnames(gas.data) = colnames(x)[c(1,5,7,9,11)]

### Create Sample used in Hamilton JOE 2003 ###
start.date = 197903
end.date = 200512
sample1 = match(start.date,x[,1]):match(end.date, x[,1])
range1  = paste(start.date,"-",end.date)

# Create lagged dependent variables
dln.rGDP.1 = lags(x[,"dln.rGDP"], 1)
dln.rGDP.2 = lags(x[,"dln.rGDP"], 2)
dln.rGDP.3 = lags(x[,"dln.rGDP"], 3)
dln.rGDP.4 = lags(x[,"dln.rGDP"], 4)
oil.data = cbind( oil.data, dln.rGDP.1, dln.rGDP.2, dln.rGDP.3, dln.rGDP.4)
gas.data = cbind( gas.data, dln.rGDP.1, dln.rGDP.2, dln.rGDP.3, dln.rGDP.4)

# Create lags of PPI_Oil and PPI_Gas
H1.PPI_Oil.1 = lags(x[,"H1.PPI_Oil"], 1)
H1.PPI_Oil.2 = lags(x[,"H1.PPI_Oil"], 2)
H1.PPI_Oil.3 = lags(x[,"H1.PPI_Oil"], 3)
H1.PPI_Oil.4 = lags(x[,"H1.PPI_Oil"], 4)
H3.PPI_Oil.1 = lags(x[,"H3.PPI_Oil"], 1)
H3.PPI_Oil.2 = lags(x[,"H3.PPI_Oil"], 2)
H3.PPI_Oil.3 = lags(x[,"H3.PPI_Oil"], 3)
H3.PPI_Oil.4 = lags(x[,"H3.PPI_Oil"], 4)

H1.PPI_Gas.1 = lags(x[,"H1.PPI_Gas"], 1)
H1.PPI_Gas.2 = lags(x[,"H1.PPI_Gas"], 2)
H1.PPI_Gas.3 = lags(x[,"H1.PPI_Gas"], 3)
H1.PPI_Gas.4 = lags(x[,"H1.PPI_Gas"], 4)
H3.PPI_Gas.1 = lags(x[,"H3.PPI_Gas"], 1)
H3.PPI_Gas.2 = lags(x[,"H3.PPI_Gas"], 2)
H3.PPI_Gas.3 = lags(x[,"H3.PPI_Gas"], 3)
H3.PPI_Gas.4 = lags(x[,"H3.PPI_Gas"], 4)

oil.data = cbind( oil.data, H1.PPI_Oil.1, H1.PPI_Oil.2, H1.PPI_Oil.3, H1.PPI_Oil.4,
                            H3.PPI_Oil.1, H3.PPI_Oil.2, H3.PPI_Oil.3, H3.PPI_Oil.4 )
gas.data = cbind( gas.data, H1.PPI_Gas.1, H1.PPI_Gas.2, H1.PPI_Gas.3, H1.PPI_Gas.4,
                            H3.PPI_Gas.1, H3.PPI_Gas.2, H3.PPI_Gas.3, H3.PPI_Gas.4 )

oil.data = data.frame(oil.data)
gas.data = data.frame(gas.data)

##### Regressions over original sample #####

### Replicate Regression 3.2 from Hamilton 2003 ###
ro3.2 = lm( dln.rGDP ~ dln.rGDP.1   + dln.rGDP.2   + dln.rGDP.3   + dln.rGDP.4   +
                       H1.PPI_Oil.1 + H1.PPI_Oil.2 + H1.PPI_Oil.3 + H1.PPI_Oil.4,
         data=oil.data, subset=sample1 )

### Replicate Regression 3.8 from Hamilton 2003 ###
ro3.8 = lm( dln.rGDP ~ dln.rGDP.1   + dln.rGDP.2   + dln.rGDP.3   + dln.rGDP.4   +
                       H3.PPI_Oil.1 + H3.PPI_Oil.2 + H3.PPI_Oil.3 + H3.PPI_Oil.4,
         data=oil.data, subset=sample1 )

### Run Regression 3.2 from Hamilton 2003 Over Gas Data ###
rg3.2 = lm( dln.rGDP ~ dln.rGDP.1   + dln.rGDP.2   + dln.rGDP.3   + dln.rGDP.4   +
                       H1.PPI_Gas.1 + H1.PPI_Gas.2 + H1.PPI_Gas.3 + H1.PPI_Gas.4,
         data=gas.data, subset=sample1 )

### Run Regression 3.8 from Hamilton 2003 Over Gas Data ###
rg3.8 = lm( dln.rGDP ~ dln.rGDP.1   + dln.rGDP.2   + dln.rGDP.3   + dln.rGDP.4   +
                       H3.PPI_Gas.1 + H3.PPI_Gas.2 + H3.PPI_Gas.3 + H3.PPI_Gas.4,
         data=gas.data, subset=sample1 )

### End Regressions over original sample ###


### Make data sets for common sample period ###
H1og.data = data.frame(na.omit(cbind(gas.data[,1:2], gas.data[,6:13], oil.data[,10:13])))
H3og.data = data.frame(na.omit(cbind(gas.data[,1:2], gas.data[,6:9],
                                     gas.data[,14:17], oil.data[,14:17])))


##### Regressions over current common full-sample #####

### Regression 3.2 from Hamilton 2003 over full data set ###
ro3.2f = lm( dln.rGDP ~ dln.rGDP.1   + dln.rGDP.2   + dln.rGDP.3   + dln.rGDP.4   +
                       H1.PPI_Oil.1 + H1.PPI_Oil.2 + H1.PPI_Oil.3 + H1.PPI_Oil.4,
         data= H1og.data )

### Regression 3.8 from Hamilton 2003 over full data set ###
ro3.8f = lm( dln.rGDP ~ dln.rGDP.1   + dln.rGDP.2   + dln.rGDP.3   + dln.rGDP.4   +
                       H3.PPI_Oil.1 + H3.PPI_Oil.2 + H3.PPI_Oil.3 + H3.PPI_Oil.4,
         data= H3og.data )

### Regression 3.2 from Hamilton 2003 Over Gas Data over full data set ###
rg3.2f = lm( dln.rGDP ~ dln.rGDP.1   + dln.rGDP.2   + dln.rGDP.3   + dln.rGDP.4   +
                       H1.PPI_Gas.1 + H1.PPI_Gas.2 + H1.PPI_Gas.3 + H1.PPI_Gas.4,
         data= H1og.data )

### Regression 3.8 from Hamilton 2003 Over Gas Data over full data set ###
rg3.8f = lm( dln.rGDP ~ dln.rGDP.1   + dln.rGDP.2   + dln.rGDP.3   + dln.rGDP.4   +
                       H3.PPI_Gas.1 + H3.PPI_Gas.2 + H3.PPI_Gas.3 + H3.PPI_Gas.4,
         data= H3og.data )


### End Regressions over current common full-sample ###

##### Restricted Regressions over current common full-sample #####

### Restricted Regressions 3.2 & 3.8 from Hamilton 2003 over full data set ###
r3.2fR = lm( dln.rGDP ~ dln.rGDP.1 + dln.rGDP.2 + dln.rGDP.3 + dln.rGDP.4, data=H1og.data )
r3.8fR = lm( dln.rGDP ~ dln.rGDP.1 + dln.rGDP.2 + dln.rGDP.3 + dln.rGDP.4, data=H3og.data )

### End Restricted Regressions over current common full-sample ###

#summary(ro3.2)
#summary(ro3.8)
#summary(rg3.2)
#summary(rg3.8)

#summary(ro3.2f)
#summary(ro3.8f)
#summary(rg3.2f)
#summary(rg3.8f)
sink(outfile)
print(anova(ro3.2f,r3.2fR))
print(anova(ro3.8f,r3.8fR))
print(anova(rg3.2f,r3.2fR))
print(anova(rg3.8f,r3.8fR))
sink()