Y<- scan("C:/CLIM/TimeS.txt")    
sink("C:/CLIM/ForecaOut.txt") ; options(digits=5)          ###
"Hohenpeissenberg, Temperature, 1781-2010"                #--#
library(TSA)                        #see CRAN software-packages
 
epsilon<- function(y,n,mc,theta,a,b){  #error terms recursively
ypr<- rep(mean(y),times=n); eps<- rep(0,times=n)  #ve. of dim n
for (t in (mc+1):n){
Suma<- theta; Sumb<- 0
for (m in 1:mc){                               
Suma<- Suma+a[m]*y[t-m]; Sumb<- Sumb+b[m]*eps[t-m]}
ypr[t]<- Suma + Sumb
eps[t]<- y[t] - ypr[t]}
return(eps)
}
                       
forec<- function(y,e,n,mc,theta,a,b,la){ #BJ forecast function
fore<- rep(mean(y),times=la)                 #vector of dim la
ep<- e ; yp<- y                              #vectors of dim n
for(j in 1:la){fore[j]<- theta
for(k in 1:mc){fore[j]<- fore[j]+a[k]*yp[n-k+1]+b[k]*ep[n-k+1]}
yp[1:(n-1)]<- yp[2:n]; yp[n]<- fore[j]
ep[1:(n-1)]<- ep[2:n]; ep[n]<- 0            #new error term = 0
}
return(fore)
}
 
#------Data---------------------------------------------------
N<- length(Y);  X<- Y                         #Y = time series 
X[1]<-0; X[2:N]<-Y[2:N]-Y[1:(N-1)]  #X=differenced series #--#
ma<- 2; mb<- 2; mc<- max(ma,mb)                  #mc maximal 6                      
c("ArOrder"=ma,"MAOrder"=mb)                  
 
# ---------- ARMA(p,q)-Model for series X --------------
xarma<- arma(X,order=c(ma,mb))
summary(xarma) 

xcoef<- xarma$coef                   #$vector of dim ma+mb+1  
#Coefficients a,b,theta 
a<- rep(0,times=6); b<- rep(0,times=6)
if (ma > 0) {for (m in 1:ma){a[m]<- xcoef[m]}}
if (mb > 0) {for (m in 1:mb){b[m]<- xcoef[ma+m]}}
theta<- xcoef[ma+mb+1]

epsil<- epsilon(X,N,mc,theta,a,b) #user function: error terms                           
 
#--------Forecasting acc. to Box&Jenkins-----------------------
la<- 10  
#user function forec: ARMA forecast function Xfore                                    
Xfore<- forec(X,epsil,N,mc,theta,a,b,la)  
"Box&Jenkins ARMA forecast vector (X series)"
Xfore; Yfore<- Xfore
 
#ARIMA forecast function Yfore (integrated series)        #--#
Yfore<- Y[N] + Xfore                                      #--#     
if(la > 1) for(j in 2:la){Yfore[j]<-Yfore[j-1]+Xfore[j]}  #--#   
"Box&J ARIMA forecast vector (Y series)"; Yfore           #--#

rm(list=objects())                                         ###