potsdTp<- read.table("C:/CLIM/PotsdT.txt",header=T)       ###
    
library(TSA)                        #see Cran-sofware-packages
attach(potsdTp)   
 
#-------Data preparation, trend removal-----------------------
"Monthly temperature means Potsdam 1893-2010"                  
mon12<- data.frame(potsdTp[,3:14])/10       #selecting jan-dec              
NYear<- length(Tyear);  M<- NYear*12
c("Number Years"=NYear,"Number Months M"=M)
detach(potsdTp)                                  

#Read the yearly trend Ytr, twelvefold as Ytre          
potsdPr<- read.table("C:/CLIM/PoT40tPr.txt",header=T)      
attach(potsdPr)                         #contains variable Ytr   
Ytre<- 1:M                          #Ytre vector of dim M  
for(m in (1:M)){j<- trunc((m-1)/12)+1; Ytre[m]<- Ytr[j]}   
##Instead of the last two lines: 
#tw<- rep(12,times=NYear); Ytre<- rep(Ytr,tw)          
YtreN1<- 9.2743       #Forecast New Year from R 4.1 (Tab. 4.2)     
 
Yobs<- as.matrix(t(mon12))               #as 12 x NYear matrix 
dim(Yobs)<- c(M,1)                       #as M-dim vector
 
Yde<- Yobs - Ytre                                  #detrending
ma<- 2; mb<- 3 ;  mc<- pmax(ma,mb)            #ARMA order <= 6  
tst<- trunc(NYear/5)*12; ts1<-tst+1      #Start for prediction  
     
#------ARMA-model for detrended monthly data------------------
#---a) Estimation-----------------------------
Ydarma<- arma(Yde,order=c(ma,mb))           
Ydcoef<- Ydarma$coef;   Yds2<- Ydarma$css
Ydfit<- Ydarma$fitted.values;    Ydres<- Ydarma$residuals
c("Start at"=tst,"ArOrder"=ma,"MAOrder"=mb,"Cond.SSQ"=Yds2) 
summary(Ydarma)

a<- rep(0,times=6); b<- rep(0,times=6)
if (ma > 0) for (m in 1:ma){a[m]<- Ydcoef[m]}
if (mb > 0) for (m in 1:mb){b[m]<- Ydcoef[ma+m]}    
theta<- Ydcoef[ma+mb+1]                           #intercept                

#---b) Prediction------------------------------  
"ARMA-prediction for the detrended series, last 12 months"   
                                       Ydfit[(M-11):M]  
YdarmaNJ<- theta
for (m in 1:mc) {YdarmaNJ<-
        YdarmaNJ+a[m]*Yde[M+1-m]+b[m]*Ydres[M+1-m]}
c("Forecast Jan_NewYear, detrended series"=YdarmaNJ)  
# (Original) series with trend
YIarma<- Ytre; YIarma[ts1:M]<- Ytre[ts1:M]+Ydfit[ts1:M] 
Yres<- Yobs - YIarma
"ARMA-prediction for the series with trend,last 12 months"  
                                       YIarma[(M-11):M] 
"Forecast Jan_NewYear, series with trend"  
c(YtreN1,YdarmaNJ,YdarmaNJ+YtreN1)                
  
#---c) Residual analysis---------------------
Yrer<- Yres[ts1:M]; msq<- mean(Yrer*Yrer)  
c("Mean Yres"=mean(Yrer),"Std Yres"=sqrt(var(Yrer)),
  "MSQ"=msq,"RootMSQ"=sqrt(msq))
racf<-acf(Yrer,lag.max=8,type="corr",plot=F)$acf    #$no Plot
"Autocorrelation of residual series"; racf[1:8]      
   
detach(potsdPr)                                        ###
rm(list=objects())                                     ###