function()
{
###
###---S+ program for Diebold, Gunther and Tay (1997)
###---Assumes database contains snp500
###---Requires module S+GARCH
###---
###---Copyright (c) 1997 F.X.Diebold, T.A.Gunther and A.S.Tay
###

#--Some utility series
	y <- snp500
	NoObs <- length(y)
	NoBins <- 20
	NoBins1 <- NoBins * 2
	Fcst <- 4133
	NoFcst <- NoObs - Fcst + 1
	l <- seq(0, 1, 0.01)	#--95% CIs for the histograms and autocorrelations
	hsehi <- rep((1 + 2 * sqrt((NoBins^2/NoFcst) * (1/NoBins) * (1 - 1/NoBins))), 101)
	hselo <- rep((1 - 2 * sqrt((NoBins^2/NoFcst) * (1/NoBins) * (1 - 1/NoBins))), 101)
	hsehi1 <- rep((1 + 2 * sqrt((NoBins1^2/NoFcst) * (1/NoBins1) * (1 - 1/NoBins1))), 101)
	hselo1 <- rep((1 - 2 * sqrt((NoBins1^2/NoFcst) * (1/NoBins1) * (1 - 1/NoBins1))), 101)
	sehi <- rep((2/sqrt(NoFcst)), 200)
	selo <- rep((-2/sqrt(NoFcst)), 200)	
#--Make forecasts and compute z's
	for(i in 1:3) {
		if(i == 1) {
			#--Naive Forecast iid N(0,1)
			z1 <- pnorm(y[Fcst:NoObs], mean = mean(y[1:(Fcst - 1)]), sd = sqrt(var(y[1:(Fcst - 1)])))
		}
		if(i == 2) {
			#--Gaussian MA(1)-GARCH(1,1)
			y.mod1 <- garch(y[1:(Fcst - 1)] ~ ma(1),  ~ garch(1, 1))
			mu <- y.mod1$coef["C", "COEF"]
			th <- ( - y.mod1$coef["MA(1)", "COEF"])
			a0 <- y.mod1$coef["A", "COEF"]
			a1 <- y.mod1$coef["ARCH(1)", "COEF"]
			b1 <- y.mod1$coef["GARCH(1)", "COEF"]
			sigma <- (a0/(1 - a1 - b1))
			temp <- 0	
			mhat <- vector("numeric", NoObs)	
			hhat <- vector("numeric", NoObs)	
			mhat[1] <- mu
			temp <- (y[1] - mu)
			hhat[1] <- a0 + b1 * sigma
			for(i in 2:NoObs) {
				mhat[i] <- mu - th * temp
				temp <- (y[i] - mhat[i])
				hhat[i] <- a0 + a1 * y[i - 1]^2 + b1 * hhat[i - 1]
			}
			sehat <- hhat^0.5	
			##---Convert to z's
			z1 <- pnorm(y[Fcst:NoObs], mhat[Fcst:NoObs], sehat[Fcst:NoObs])
		}
		if(i == 3) {
			#--t-MA(1)-GARCH(1,1)
			y.mod1 <- garch(y[1:(Fcst - 1)] ~ ma(1),  ~ garch(1, 1), cond.dist = "t")
			mu <- y.mod1$coef["C", "COEF"]
			th <- ( - y.mod1$coef["MA(1)", "COEF"])
			vv <- y.mod1$cond.dist$dist.par
			a0 <- y.mod1$coef["A", "COEF"]
			a1 <- y.mod1$coef["ARCH(1)", "COEF"]
			b1 <- y.mod1$coef["GARCH(1)", "COEF"]
			sigma <- (a0/(1 - a1 - b1))
			temp <- 0
			mhat <- vector("numeric", NoObs)
			hhat <- vector("numeric", NoObs)
			mhat[1] <- mu
			temp <- (y[1] - mu)
			hhat[1] <- a0 + b1 * sigma
			for(i in 2:NoObs) {
				mhat[i] <- mu - th * temp
				temp <- (y[i] - mhat[i])
				hhat[i] <- a0 + a1 * y[i - 1]^2 + b1 * hhat[i - 1]
			}
			sehat <- hhat^0.5
			temp <- (y[Fcst:NoObs] - mhat[Fcst:NoObs])/sehat[Fcst:NoObs]
			##--Convert to z's
			z1 <- pt(temp, vv)
		}
		##--Computing acfs
		z2 <- (z1 - mean(z1))^2
		z3 <- (z1 - mean(z1))^3
		z4 <- (z1 - mean(z1))^4
		z1.acf <- acf(z1, lag.max = 200, type = "correlation", plot = F)
		z2.acf <- acf(z2, lag.max = 200, type = "correlation", plot = F)
		z3.acf <- acf(z3, lag.max = 200, type = "correlation", plot = F)
		z4.acf <- acf(z4, lag.max = 200, type = "correlation", plot = F)
		z1.acf <- z1.acf$acf[2:201, 1, 1]
		z2.acf <- z2.acf$acf[2:201, 1, 1]
		z3.acf <- z3.acf$acf[2:201, 1, 1]
		z4.acf <- z4.acf$acf[2:201, 1, 1]

		#--plotting histograms
		win.graph(width = 170, height = 90, density = 10, units.per.inch = 25.4)
		if(i == 1) {
			par(fig = c(0, 0.5, 0, 1), mar = c(2, 2, 2, 2))
			par(cex = 0.7, mex = 0.7)
			hist(z1, NoBins, col = 0, ylim = c(0, 2), ylab = "", xlab = "", axes = T, xaxs = "s", yaxs = "s", 
				probability = T)
			ylo <- par()$yaxp[1]
			yhi <- par()$yaxp[2]
			box()
			par(new = T)
			plot(l, hsehi, type = "l", lty = 2, ylim = c(0, 2), ylab = "", xlab = "", axes = T, xaxs = "s", yaxs
				 = "s")
			par(new = T)
			plot(l, hselo, type = "l", lty = 2, ylim = c(ylo, yhi), ylab = "", xlab = "", axes = T, xaxs = "s", 
				yaxs = "s")
			par(fig = c(0.5, 1, 0, 1), mar = c(2, 2, 2, 2))
			par(cex = 0.7, mex = 0.7)
			hist(z1, 40, col = 0, ylim = c(0, 2), ylab = "", xlab = "", axes = T, xaxs = "s", yaxs = "s", 
				probability = T)
			ylo <- par()$yaxp[1]
			yhi <- par()$yaxp[2]
			box()
			par(new = T)
			plot(l, hsehi1, type = "l", lty = 2, ylim = c(0, 2), ylab = "", xlab = "", axes = T, xaxs = "s", 
				yaxs = "s")
			par(new = T)
			plot(l, hselo1, type = "l", lty = 2, ylim = c(ylo, yhi), ylab = "", xlab = "", axes = T, xaxs = "s", 
				yaxs = "s")
		}
		if(i == 2 || i == 3) {
			par(fig = c(0, 1, 0, 1), mar = c(3, 3, 3, 3))
			par(cex = 0.7, mex = 0.7)
			hist(z1, NoBins, col = 0, ylim = c(0, 2), ylab = "", xlab = "", axes = T, xaxs = "s", yaxs = "s", 
				probability = T)
			ylo <- par()$yaxp[1]
			yhi <- par()$yaxp[2]
			box()
			par(new = T)
			plot(l, hsehi, type = "l", lty = 2, ylim = c(0, 2), ylab = "", xlab = "", axes = T, xaxs = "s", yaxs
				 = "s")
			par(new = T)
			plot(l, hselo, type = "l", lty = 2, ylim = c(ylo, yhi), ylab = "", xlab = "", axes = T, xaxs = "s", 
				yaxs = "s")
		}

		##--plotting acfs
		win.graph(width = 180, height = 130, density = 10, units.per.inch = 25.4)
		par(cex = 0.7, mex = 0.7, fig = c(0, 0.5, 0.5, 1), mar = c(3, 3, 3, 3))
		plot(z1.acf, type = "h", ylim = c(-0.1, 0.5), xlab = "", ylab = "")
		par(new = T)
		plot(sehi, yaxs = "d", xaxs = "d", type = "l", lty = 2, xlab = "", ylab = "")
		par(new = T)
		plot(selo, yaxs = "d", xaxs = "d", type = "l", lty = 2, xlab = "", ylab = "")
		mtext("(a)", side = 1, line = 2)
		par(fig = c(0.5, 1, 0.5, 1), mar = c(3, 3, 3, 3))
		plot(z2.acf, type = "h", ylim = c(-0.1, 0.5), xlab = "", ylab = "")
		par(new = T)
		plot(sehi, yaxs = "d", xaxs = "d", type = "l", lty = 2, xlab = "", ylab = "")
		par(new = T)
		plot(selo, yaxs = "d", xaxs = "d", type = "l", lty = 2, xlab = "", ylab = "")
		mtext("(b)", side = 1, line = 2)
		par(fig = c(0, 0.5, 0, 0.5), mar = c(3, 3, 3, 3))
		plot(z3.acf, type = "h", ylim = c(-0.1, 0.5), xlab = "", ylab = "")
		par(new = T)
		plot(sehi, yaxs = "d", xaxs = "d", type = "l", lty = 2, xlab = "", ylab = "")
		par(new = T)
		plot(selo, yaxs = "d", xaxs = "d", type = "l", lty = 2, xlab = "", ylab = "")
		mtext("(c)", side = 1, line = 2)
		par(fig = c(0.5, 1, 0, 0.5), mar = c(3, 3, 3, 3))
		plot(z4.acf, type = "h", ylim = c(-0.1, 0.5), xlab = "", ylab = "")
		par(new = T)
		plot(sehi, yaxs = "d", xaxs = "d", type = "l", lty = 2, xlab = "", ylab = "")
		par(new = T)
		plot(selo, yaxs = "d", xaxs = "d", type = "l", lty = 2, xlab = "", ylab = "")
		mtext("(d)", side = 1, line = 2)
	}
	z1
}