#Hohenpeissenberg, Temperatur, Jahr/Winter, mit Polynom Ausgleich
# R: text( ) nicht ausserhalb des Plot-Rahmens
#Autokorr und Partial Autocorr
#Periodogramm und Spektraldichte

hohtp<- read.table("Hhohen.dat",header=T)         #evtl: Hhohen.txt
attach(hohtp)
 
postscript(file="HPja0.ps",height=18,width=10,horizontal=F)
par(mfrow=c(2,1))

Y<- Tjahr/100
length(Y); mean(Y); sqrt(var(Y))
plot(Jahr,Y,type="l",lty=1,xlim=c(1780,2005),ylim=c(4.0,9.0),
 xlab="Jahr",ylab="Mittlere Jahrestemperatur")
title(main="Hohenpeissenberg, Temperatur Jahr 1781-2005",cex=0.85)
abline(h=c(4,5,8),lty=2)
abline(h=6.333,lty=3)
text(1880,8.9,"Jahres-Mittelwerte (mit Ausgleichs-Polynom 4. Grades)",cex=0.8) 
text(2025,6.33,"6.33",cex=0.7)

Ja<- (Jahr-1780)/100
J2<- Ja*Ja
J3<- J2*Ja
J4<- J3*Ja
tppol<- lm(Y~Ja+J2+J3+J4)
#summary(tppol)
lines(Jahr,predict(tppol),lty=2)

#Winter
Twi<- (dzvj+jan+feb)/30

plot(Jahr,Twi,type="l",lty=1,xlim=c(1780,2005),ylim=c(-6.5,3.5),
   xlab="Jahr",ylab="Mittlere Wintertemperatur")
title(main="Hohenpeissenberg, Temperatur Winter 1781-2005",cex=0.85)
abline(h=c(-6,-4,,2),lty=2)
abline(h=-1.368,lty=3)
text(1880,3.1,"Winter-Mittelwerte (mit Ausgleichs-Polynom 4.Grades)",cex=0.8)
text(2025,-1.37,"-1.37",cex=0.7)

Ja<- Jahr -1780
Y<-Twi
length(Y); mean(Y); sqrt(var(Y)) 
J2<- Ja*Ja
J3<- J2*Ja
J4<- J3*Ja
tppol<- lm(Y~Ja+J2+J3+J4)
#summary(tppol)
lines(Jahr,predict(tppol),lty=2)

dev.off()

#Residuen vom Polynom-Trend (Trendbereinigung der Winterreihe)
Yres<- Y - predict(tppol)
c(mean(Yres),sqrt(var(Yres)))

postscript(file="HPja1.ps",height=18,width=10,horizontal=F)
par(mfrow=c(2,1))

#Autocorr und Partial Autocorr fuer Residuen-Zeitreihe
acf(Yres,lag.max=20,type="corr")
title(main="Autocorr und Partial Autocorr: TpWinter-Residuen",cex=0.7)
acf(Yres,lag.max=20,type="partial")

dev.off()

postscript(file="HPja2.ps",height=18,width=10,horizontal=F)
par(mfrow=c(2,1))

#Periodogramm und Spektraldichte fuer Residuen-Zeitreihe
hwsp1<- spectrum(Yres,detrend=F,demean=T,plot=T,ylab="Periodogramm")
title(main="Periodogramm+Spektraldichte: TpWinter-Residuen",cex=0.7)
hwsp2<- spectrum(Yres,spans=21,detrend=F,demean=T,plot=T,            
   ylab="Spektraldichte (span=21)")     
                 
dev.off()

rm(Y,Ja,J2,J3,J4,tppol)
rm(Twi)