bremenTp<- read.table("C:/CLIM/BremenT.txt",header=T)         ###  
 
attach(bremenTp)
quot<- "Bremen, Temperature, 1890-2010"; quot
 
library(TSA)                         #see Cran-Software-Packages
#---------------------------------------------------------------
armat<- function(Y,n,tst,ma,mb){   #forecast regression approach
#parmat vect of dim n+1,components 1,..,tst filled with mean val
parmat<- 1:(n+1); parmat[1:tst]<- rep(mean(Y[1:tst]),times=tst)
 
for(t in tst:n) {
art<- arma(Y[1:t],order=c(ma,mb))      
coef<- art$coef;   resi<- art$residuals
parma<- coef[ma+mb+1]                           #intercept theta
a<- rep(0,times=12); b<- a; mc<- pmax(ma,mb) 
if (ma > 0) for (m in 1:ma){a[m]<- coef[m]}         #alpha-coeff.
if (mb > 0) for (m in 1:mb){b[m]<- coef[ma+m]}       #beta-coeff.
for (m in 1:mc){parma<- parma + a[m]*Y[t+1-m] + b[m]*resi[t+1-m]}
parmat[t+1]<- parma
} 
return(parmat)  }                #return prediction vector parmat

#---------------Data preparation, Differencing-------------------
mon12<- data.frame(bremenTp[,3:14])/10;           #select jan-dec                  
Yje<- rowMeans(mon12)          #more precise than Yje<- Tyear/100      
N<- length(Year);                 Dy<- Yje                                           
Dy[1]<- 0; Dy[2:N]<- Yje[2:N]-Yje[1:(N-1)] #Dy differenced series 

ma<- 3; mb<- 1; tst<- trunc(N/5);  ts1<- tst+1          
c("ArOrder"=ma,"MAOrder"=mb,"Start"=tst,"Ncut"=N-tst,
                "StdevDYcut"=sqrt(var(Dy[ts1:N])))                  
              
# ---------- ARIMA(p,q)-model and ARIMA-prediction----- ---------
#  a) Differenced series Dy---------------------
darma<- arma(Dy,order=c(ma,mb));  "Results for whole diff series"
summary(darma)                                   #for output only
"Forecast regression approach, differenced series"
parmat<- armat(Dy,N,tst,ma,mb)                   #vector i=1:(N+1) 
pred<- parmat[ts1:N]; dres<- Dy[ts1:N]-pred  #vector i=1:(N-(tst))
msq<- mean(dres*dres)
c("MeanDres"=mean(dres),"StdDres"=sqrt(var(dres)),
  "MSQ (differenced)"=msq,"RootMSQ"=sqrt(msq))         
"Auto-correlations of residuals, differenced series"
acf(dres,lag.max=8,type="corr",plot=F)
 
#  b) Integrated series YIarma----------------------
YIarma<- Yje; YIarma[2:N]<- Yje[1:(N-1)]+parmat[2:N]  #vect i=1:N
"Observations and ARIMA-predictions for last decade"  
Yje[(N-9):N]; YIarma[(N-9):(N)]
c("Forecast NewYear"= Yje[N],parmat[N+1],Yje[N]+parmat[N+1])    

#--Supplement---------------------------------------
"Test: ARIMA-residuals = ARMA-residuals of differenced series" 
"Test: same MSQ and auto-corr values as above"
YIres<- Yje - YIarma                   #YIres[ts1..] = dres[1...] 
YIrer<- YIres[ts1:N];  misq<- mean(YIrer*YIrer)      
c("MSQ (integrated)"=misq,"RootMSQ"=sqrt(misq)) 
"Auto-correlations of ARIMA-residuals"
acf(YIrer,lag.max=8,type="corr",plot=F) 

detach(bremenTp)                                              ###                                             
rm(list=objects())                                            ###