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)
}
                       
monca<- function(y,e,n,mc,theta,a,b,la,se){      #MC simulation
monc<- rep(mean(y),times=la)
ep<- e ; yp<- y; eps<- rnorm(la)      #la N(0,1) random numbers
for(j in 1:la){monc[j]<- theta+eps[j]*se
for(k in 1:mc){monc[j]<- monc[j]+a[k]*yp[n-k+1]+b[k]*ep[n-k+1]}
yp[1:(n-1)]<- yp[2:n]; yp[n]<- monc[j]
ep[1:(n-1)]<- ep[2:n]; ep[n]<- eps[j]*se        #new error term
}
return(monc)
}

#------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                           
se<- sqrt(var(epsil[(mc+1):N]))    #error terms from mc+1 onw.
c("Mean Epsilon"=mean(epsil[(mc+1):N]),  "StdDev Epsilon"=se)
 
#--------Forecasting acc. to Monte Carlo-----------------------
la<- 10; MC<- 20000; c("Monte Carlo Repetitions"=MC)  
Xmonca<- 1:(MC*la); dim(Xmonca)<-c(MC,la)  #MCxla matrix Xmonca
#user function monca: 1 Monte-Carlo repetition
for (m in 1:MC){ 
montc<- monca(X,epsil,N,mc,theta,a,b,la,se) 
Xmonca[m,]<- montc }
Ymonca<- Xmonca

#MonteCarlo simulations Ymonca (integrated series)         #--#
Ymonca<- Y[N] + Xmonca                                     #--#
if(la>1) for(j in (2:la))                                  #--#
                {Ymonca[,j]<- Ymonca[,(j-1)]+Xmonca[,j]}   #--#
  
meYmonca<- colMeans(Ymonca)
"Monte Carlo mean vector (Y series)";  meYmonca
quan<- 1:(4*la); dim(quan)<- c(4,la)          #4xla matrix quan
alph<- c(0.10,0.20,0.80,0.90)
for(j in 1:la){
quan[,j] <- quantile(Ymonca[,j],alph) }    #empirical quantiles

"Monte Carlo quantile vectors, levels 0.10, 0.20, 0.80, 0.90"
quan[,] 

rm(list=objects())                                         ###