postscript(file="C:/CLIM/Intfunc.ps",height=6,width=16,horiz=F)

#----------------------------------------------------------------
gauss<- function(t,x,b){                            #Gauss kernel
c<- sqrt(2*pi)
fun<- (1/c)*exp(-(t-x)^2/(2*b*b)) 
return(fun)
}
rectan<- function(t,x,b){                     #Rectangular kernel
fun<- 0
if({t-b <= x} & {x<= t+b}) fun<- 1/2
return(fun)
}
 
#-----------------------------------------------------------------
quot<- "Aachen 2004-2010. Daily Precipitation"; quot

xx<- scan("C:/CLIM/eventAa.txt")
n<- xx[1];  tx<- xx[2:(n+1)];  tb<- 2555   
c("n"=n,"Right end"=tb,"lambda"=n/tb,"log lambda"=log(n/tb))

"Nonparametric kernel estimation"
ndelt<-200; delt<- tb/ndelt; bh<- 40                 #bh bandwidth                    
tt<- delt*(1:ndelt); int<- 1:ndelt           #vectors of dim ndelt

for(j in 1:ndelt){ kern<- 0  
for(i in 1:n){                     #choose: gauss( ) or rectan( ):
kern<- kern + gauss(tt[j],tx[i],bh)}   
int[j]<- (1/bh)*kern}                            

d<- bh; red<-tt>d & tt<tb-d       #logical vector, margin width d   
ttr<-tt[red]; intr<-int[red] 
c("ndelt"=ndelt,"delta"=delt,"bandwith"=bh,
                      "mean int"=mean(intr))
 
plot(ttr,intr,type="l",lty=1,xlab="Days (Occurrence times +)",
                       ylab="Nonparametric Intensity Function")
title(main=quot)  
abline(h=n/tb,lty=2)
text(tb+10,n/tb,"l",font=5)                        #Greek lambda
abline(v=365*(1:7),lty=3)
text(tx,min(intr),"+",cex=0.8)            #mark occurrence times

dev.off()
rm(list=objects())                                           ###