n<-340

# Pseudo Simulate an ARCH(1) model

alpha0<-500
garcherr<-function(n){rchisq(n,1)-1}
garchsq<-arima.sim(list(ar=c(.5)),n=n,rand.gen=garcherr)+alpha0
garchp<-garchsq^.5
par(mfrow=c(2,1))
tsplot(garchp)
title(main="Pseudo Simulated GARCH Process")
tsplot(garchsq)
title(main="Squared GARCH Process")

par(mfrow=c(2,2))
hist(garchp)
hist(garchsq)
acf(garchsq)
acf(garchsq,type="partial")


# Fit the parameters
garchapprox<-ar.gm(garchsq,order=10)

# Quick check on approximate approach

mtheory<-garchapprox$rmu/(1-garchapprox$ar)
mdata<-mean(garchsq)

# Let's take a look at the ACF 
# How could we get an estimate of the PACF?

par(mfrow=c(2,1))
acf(garchsq)
tsplot(garchapprox$chat[1,]/garchapprox$chat[1,1], type="p")

title(main="Estimator of ACF Based on Robust AR fit")

# Compute one-step ahead predictions for garchsq - based on L2
gforecast<-arima.filt(garchsq,list(ar=garchapprox$ar,sigma2=garchapprox$sd),
	xreg=(garchsq*0+1),reg.coef=garchapprox$rmu)
	
par(mfrow=c(1,1))	
tsplot(gforecast$filt,gforecast$pred)

# Compute and plot standardized residuals 

par(mfrow=c(2,2))

gresid<-(gforecast$filt-gforecast$pred)/(gforecast$pred^.5)
tsplot(gresid)
plot(gforecast$filt,gforecast$pred)
hist(gresid)
boxplot(gresid)

# NOte there is no problem in generalizing the above to a GARCH(1,0)