library(TSA)                        #see Cran-software-packages
Htpres<- read.table("C:/CLIM/HoT22Res.txt",header=T)  
attach(Htpres) 
quot<- "Hohenpberg, Temp 1781-2010, Residuals from trend"; quot           
N<-  length(Year); sde<- sqrt(var(Y))                          
c("Mean Y"=mean(Y),"Stdev Y"=sde,"Number Years"=N)             
      
#---------R function GARCH-------------------------------------
ma<- 3;  mb<- 1 ;   mc<- pmax(ma,mb)+1        #Order ma,mb <= 4  
c("ma"=ma, "mb"=mb, "mc"=mc)                           
ord<- c(mb,ma)                            #reversed order input
zgarch<- garch(Y,order=ord,maxiter=40)             
summary(zgarch)                          #a.o. diagnostic tests                        
         
 #---------GARCH-prediction, sigma(t) estimation---------------
"Spre GARCH-prediction for sigma(t)" 
Spre<- zgarch$fitted.values                      #$vector dim N
"GARCH(p,q)-prediction, last ten values"; Spre[(N-9):N]  
#Shat emp. estimator for sigma(t), moving blocks of length ka
ka<- 5; Shat<- rep(1,times=N)                
for(t in ka:N) {Shat[t]<- sqrt(var(Y[(t-ka+1):t]))} 
"Sigma(t)-estimation, last ten values"; Shat[(N-9):N]  

#----------GARCH-residuals epsilon  eps(t)----------------------
res<- Y/Spre; eps<- res[mc:N]
c("Mean epsilon"=mean(eps),"Stdev epsilon"=sqrt(var(eps)))
racf<- acf(eps,lag.max=8,type="corr",plot=F)
"Auto-correlations of GARCH-residuals epsilon"; racf$acf[1:8] #$
  
#--------Plot---------------------------------------------------
postscript(file="C:/CLIM/GARCHmod.ps",height=6,width=16,horiz=F) 
cylim<-c(-1.5,1.5); ytext<- "Temperature [C]"                  
plot(Year,Y,type="l",lty=1,xlim=c(1780,2010),              
                    ylim=cylim, xlab="Year",ylab=ytext,cex=1.3)
title(main=quot) 
lines(Year[mc:N],Spre[mc:N],type="l",lty=1)
lines(Year[mc:N],Shat[mc:N],type="l",lty=2)
abline(h=c(-1,0,1),lty =3)
 
dev.off()                                                              
#---- Supplement ---------------------------
 
garchpr<- function(y,n,a,ma,b,mb,sde){         #GARCH prediction
ys<- rep(sde,times=n)                           #vector of dim n                  
mc<- pmax(ma,mb)+1                              
for (t in mc:n){
suma<- a[1]; sumb<- 0                        #a[1] constant term
for (m in 1:(mc-1)){                               
suma<- suma + a[m+1]*y[t-m]^2
sumb<- sumb + b[m]*ys[t-m]^2}
ys[t]<- sqrt(suma+sumb)}
return(ys)                          #return prediction vector ys
}

a<- rep(0,times=5); b<- rep(0,times=4)
for (m in 1:(ma+1))  {a[m]<- zgarch$coef[m]}    #$a[1] constant    
if (mb > 0) {for (m in 1:mb)  {b[m]<-zgarch$coef[m+ma+1]}}   #$ 

"Calculation of GARCH-prediction per user function garchpr" 
spre<- garchpr(Y,N,a,ma,b,mb,sde)          
"GARCH(p,q)-prediction, last ten values, user function"
spre[(N-9):N]                                 #spre[ ] = Spre[ ]
  
detach(Htpres)                                              ###       
rm(list=objects())                                          ###