# gdp2.r #

library(quadprog)		# You may need to install the package #
library(tseries)		# You may need to install the package #

gdpdata <- read.table("gdpdata.txt")
gdpdata <- as.matrix(gdpdata)
n=nrow(gdpdata)
year=as.matrix(gdpdata[1:n,1])
quarter=as.matrix(gdpdata[1:n,2])
t=year+(quarter-1)/4

gdp=as.matrix(gdpdata[1:n,3])
gdpg=as.matrix(gdpdata[1:n,4])
t10=as.matrix(gdpdata[1:n,5])
t3=as.matrix(gdpdata[1:n,6])
aaa=as.matrix(gdpdata[1:n,7])
baa=as.matrix(gdpdata[1:n,8])
hs=as.matrix(gdpdata[1:n,9])
bp=as.matrix(gdpdata[1:n,10])
spread=t10-t3
dspread=baa-aaa

# Create Data Matrix Using 12 Initial Conditions @
kk=12			# Number of initial conditions #
nn=n-kk			# Number of data points less number of initial conditions equals number of observations #
y=as.matrix(gdpg[(1+kk):n])	# dependent variable #
x=matrix(1,nn+1,1)	# Regressors, first column (ones), one more observation than dependent variable (for forecast) #
for (j in 1:kk)
{ x=cbind(x,gdpg[(1+kk-j):(n-j+1)]) } 	# X matrix columns, lags of y #


# Model Combination #
kn=kk+1				# Number of models = AR(0) through AR(kk) #
yf=matrix(0,kn,1)		# vector of forecasts (empty for now) #
ee=matrix(0,nn,kn)		# matrix of prediction errors (empty for now) #

for (k in 1:kn)
{
  xk=x[1:nn,1:k]
  xf=x[nn+1,1:k]
  xxi=solve(t(xk)%*%xk)
  beta=xxi%*%(t(xk)%*%y)
  e=y-xk%*%beta
  h=rowSums((xk%*%xxi)*xk)
  eh=e/(1-h)
  yf[k]=xf%*%beta  
  ee[,k]=eh	
}

Dmat=(t(ee)%*%ee)/nn
dvec=matrix(0,kn,1)
Amat=t(rbind(matrix(1,1,kn),diag(kn)))
bvec=rbind(1,matrix(0,kn,1))
QP <- solve.QP(Dmat,dvec,Amat,bvec,bvec)
w <- QP$solution
w <- as.matrix(w)
e=ee%*%w
cv=t(w)%*%Dmat%*%w

yff=t(yf)%*%w
print("Models, Weights")
print(cbind(seq(0,kk),w))
print("CV, Combination Forecast")
print(cbind(cv,yff))


# Leading Indicator Forecasting MOdel #
xs=cbind(x[1:(nn+1),1:3],spread[(kk):n],dspread[(kk):n],hs[(kk):n],bp[(kk):n]) # Leading Indicators #

xn=ncol(xs)
s=seq(1,xn)		# column indicators for xs #
			# We now create a matrix, where each row indicates which elements of xs to include in a model #
			# Each row is a model #
			# The number of columns is the same as xs #
models1=c(		
1,1,0,0,0,0,0,
1,1,0,1,0,0,0,
1,1,0,0,1,0,0,
1,1,0,0,0,1,0,
1,1,0,0,0,0,1,
1,1,0,1,1,0,0,
1,1,0,1,0,1,0,
1,1,0,1,0,0,1,
1,1,0,0,1,1,0,
1,1,0,0,1,0,1,
1,1,0,0,0,1,1,
1,1,0,1,1,1,0,
1,1,0,1,1,0,1,
1,1,0,1,0,1,1,
1,1,0,0,1,1,1,
1,1,0,1,1,1,1,
1,1,1,0,0,0,0,
1,1,1,1,0,0,0,
1,1,1,0,1,0,0,
1,1,1,0,0,1,0,
1,1,1,0,0,0,1,
1,1,1,1,1,0,0,
1,1,1,1,0,1,0,
1,1,1,1,0,0,1,
1,1,1,0,1,1,0,
1,1,1,0,1,0,1,
1,1,1,0,0,1,1,
1,1,1,1,1,1,0,
1,1,1,1,1,0,1,
1,1,1,1,0,1,1,
1,1,1,0,1,1,1,
1,1,1,1,1,1,1
)

jj=length(models1)/xn			# number of models #
# matrix of variables to select #
models=matrix(models1,nrow=jj,ncol=xn,byrow=1)

yf=matrix(0,jj,1)			# vector of forecasts (empty for now) #
ee=matrix(0,nn,jj)		# matrix of prediction errors (empty for now) #
				
for (j in 1:jj){
  ji=s[models[j,]==1]
  xj=xs[,ji]
  xk=xj[1:nn,]
  xf=xj[nn+1,]
  xxi=solve(t(xk)%*%xk)
  beta=xxi%*%(t(xk)%*%y)
  e=y-xk%*%beta
  h=rowSums((xk%*%xxi)*xk)
  eh=e/(1-h)
  ee[,j]=eh
  yf[j]=xf%*%beta  
}

Dmat=(t(ee)%*%ee)/nn
dvec=matrix(0,jj,1)
Amat=t(rbind(matrix(1,1,jj),diag(jj)))
bvec=rbind(1,matrix(0,jj,1))
QP <- solve.QP(Dmat,dvec,Amat,bvec,bvec)
w <- QP$solution
w <- as.matrix(w)
e=ee%*%w
yff=t(yf)%*%w
cv=t(w)%*%Dmat%*%w
sig=(t(e)%*%e)/nn

print("Model, Weights")
print(cbind(models,w))
print("CV, Combination Forecast")
print(cbind(cv,yff))

st=c(.8,.2,.7)
x.arch <- garch(e,order=c(1,1),control=garch.control(start=st))
print(summary(x.arch))
archc=coef(x.arch)
sd=predict(x.arch)
like=logLik(x.arch)

var <- as.matrix(sd[,1]^2)
varf=archc[1]+archc[2]*(e[nn]^2)+archc[3]*var[nn,1]
print("Forecast Variance, Standard Deviation")
print(cbind(varf,sqrt(varf)))
print("Unconditional Variance, Standard Deviation")
print(cbind(sig,sqrt(sig)))
print("Log Likelihood")
print(like)

t1 <- as.matrix(t[(kk+1):n])
e2=e^2
plot(t1,e,main="Leave-One-Out Prediction Residuals",type="l",xlab="",ylab="")
windows()
plot(t1,e2,main="Squared Prediction Residuals",type="l",xlab="",ylab="")
windows()
plot(t1,var,main="Estimated Variance",type="l",xlab="",ylab="")