# ----------------- # Question 1 data <- read.table("../Data/Compstat/investment.txt",header=TRUE) attach(data) library("MASS") y<- Invest/(CPI*10) time<- 1:15 GNP1<-data$GNP/(CPI*10) X<- cbind(rep(1,15),time,GNP1,Interest,Inflation) data1<-data.frame(Invest=y,time, GNP=GNP1,Interest,Inflation) plot(data1) reg<- lm(y~time+GNP1+Interest+Inflation) # Least-squares estimation # ----------------- # Question 2 # Gibbs sampler gibbs_reg<- function(y,X,beta0,B0,alpha0,delta0,B,N){ n<-dim(X)[1] p<-dim(X)[2] B0inv<-solve(B0) alpha1=alpha0+n beta<-matrix(0,nrow=B+N,ncol=p) sigma2<-vector(length=N,mode="numeric") # Initialization sigma20<-1/rgamma(1, shape=alpha0/2, rate = delta0/2) for(t in 1:(B+N)){ print(t) B1<-solve(sigma20^{-1}* t(X)%*%X + B0inv) beta1<- B1 %*% (sigma20^{-1}*t(X)%*%y + B0inv%*%beta0) delta1<- delta0+t(y-X%*%beta1)%*%(y-X%*%beta1) beta[t,]<- mvrnorm(n=1,mu=beta1,Sigma=B1) sigma2[t]<- 1/rgamma(1, shape=alpha1/2, rate = delta1/2) sigma20<-sigma2[t] } return(list(beta=beta[(B+1):(B+N),],sigma2=sigma2[(B+1):(B+N)])) } # ----------------- # Question 3 p<-dim(X)[2] V<-1000 E<-0.01 alpha0<-2*(E^2/V+2) delta0<-2*E*(E^2/V+1) B0<-10*diag(p) beta0<-rep(0,p) B<- 1000 N<-5000 sim<-gibbs_reg(y,X,beta0,B0,alpha0,delta0,B,N) # Diagnostics par(mfrow=c(1,2)) plot(sim$sigma2,type="l") acf(sim$sigma2,lag.max=200) par(mfrow=c(1,1)) par(mfrow=c(1,2)) plot(sim$beta[,1],type="l") acf(sim$beta[,1],lag.max=200) par(mfrow=c(1,1)) hist(sim$sigma2) # ----------------- # Question 4 predict_reg<- function(x0,sim){ N<- length(sim$sigma2) Ey0<- vector(length=N,mode="numeric") y0<-Ey0 for(t in 1:N){ Ey0[t]<- t(x0)%*%sim$beta[t,] y0[t] <- Ey0[t] + sqrt(sim$sigma2[t])*rnorm(1) } return(list(Ey0=Ey0,y0=y0)) } x0<-c(1,16,1.6,12,7) pred<-predict_reg(x0,sim) par(mfrow=c(1,2)) hist(pred$Ey0) hist(pred$y0) par(mfrow=c(1,1)) boxplot(pred$Ey0,pred$y0)