setwd("C:/Dokumente und Einstellungen/zeugner/Eigene Dateien/IHS/bayes")  #set working directory, 
#pay attention to / instead of the common backslash
#install packages MCMCpack and lmtest via the menu
library(MCMCpack)  #load the package for Gibbssampling
library(lmtest)  # load the package for Dirty Frequentist tests
library(xtable) #for generating latex output tables
gasoline<-read.table("gasoline.txt",header=T) #loading the data, header means first row contains variable names, T means True
attach(gasoline) #attaches the data frame to be able to work with variable names directly 
summary(gasoline)
G.pop <- Gascon/Pop  #defines a new object (variable) G.pop which is the Gasconsumption divided by population
ls() #to see which objects are currently stored in your workspace
rm(G.pop) #deletes the G.pop object, rm(list=ls()) deletes everything in the workspace


########
#matrices, vectors
#######
t<-c(1:36)  #time trend, c defines a row vector with values 1, 2, 3, .....,36
a<-rep(4,3) #defines a vector a=(4,4,4)
c<-matrix(c(3,2,1,1,2,3),nrow=2,byrow=T)  #defines a 2 times 3 matrix
b<-c[1,1] #selects the first element of the matrix
b<-c[,1]  #all rows, first column
b<-c[1,] #first row, all columns
a<-cbind(c[,1],c[,2])  #defines a new object consisting of the first column vector of c and the second one
a<-rbind(c[1,],c[2,])  #rowbind instead of columnbind

#in general, type ?"command" to get help for your command, for example ?plot, also in the
#tool menu, look at html help, then click on the package you need to get a list of all its commands with examples

######################################
#OLS Regression as in Chapter 2, p 12
#######################################
plot(Gascon,Pgas)
#ols regression with constant, see ?lm how to regress without a constant
model1.ols<-lm(log(G.pop) ~ log(Y) + log(Pgas) + log(Pncar) + log(Pucar))  
summary(model1.ols)
names(model1.ols) # to get the names of the objects stored in model1.ols

coeff<-model1.ols$coefficients #saves the previously estimated coefficients in a vector
coeff<-model1.ols[1] #equivalently

model.output<-xtable(model1.ols)  #for latex user gives you the regression table in latex format
print(model.output)

model2.ols<-lm(log(G.pop) ~ log(Y/Pop) + log(Pgas) + log(Pncar) + log(Pucar) + Year)
summary(model2.ols)
dwtest(log(G.pop) ~ log(Y/Pop) + log(Pgas) + log(Pncar) + log(Pucar) + Year)

###########################################
#Koop section 4 & MCMC
###############################################
hprice<-read.table("HPRICE.txt",header=T)
attach(hprice)

#koop informative prior
koop2<-MCMCregress(formula=percentprice ~ lot.siz + nrbed + nrbath  + nrstor, data = hprice, burnin = 1000, mcmc= 10000,
verbose = 0, seed = NA, beta.start = NA, b0= c(0,10,5000,10000,10000),B0=c(10000^-2,5^-2,2500^-2,5000^-2,5000^-2),c0=5,d0=5*((4e-08)^-1),marginal.likelihood ="Chib95")


#b0 is the prior vector like c(0.1,0,1)
#B0 is the inverse V from koop, the precision prior, you enter the diagonal elements
# for example c(2^-2,3^-2,1^-2)


summary(koop2)
plot(koop2)
HPDinterval(koop2)
geweke.diag(koop2)  #geweke cd
geweke.plot(koop2)

?BayesFactor





















Done by DER_FLIEGA