###############################################################
### evs_assignments1.R
### Nonlinear data analysis and modeling in sciences, SS07
### Extreme Value Statistics
### Henning Rust                                     28.06.2007


### read the data file
file <- "evs_assignments1.dat"                                 ## define file name
dat <- read.table(file,col.names=c("day","month","year","Q"))  ## read data
dat.ts <- ts(dat$Q,start=c(305,1850),frequency=365)            ## convert to time series


### plot the time series
plot(dat.ts,ylab="discharge [m^3/s]",main="River Discharge")   ## plot ts


### histogram
hist(dat.ts,xlab="discharge [m^3/s]",main="River Discharge",sub="60 equally spaced breaks",xlim=c(0,2000),breaks=60,prob=FALSE)


### kernel estimates of the probability density
plot(density(dat$Q,bw=20,kernel="gaussian"),
     xlab="discharge [m^3/s]",main="River Discharge",sub="Kernel: Gaussian, sd=20",xlim=c(0,2000))


### estimate parameters of a log-normal distribution by methods of moments
log.dat <- log(dat.ts)                     ## log-transform
log.mu <- mean(log.dat)                    ## estimate mean of transformed data 
log.sd <- sd(log.dat)                      ## estimate standard deviation of transformed data
x <- seq(0,2500,length=500)                ## generate a sequence of quantiles
hist(dat$Q,xlab="discharge [m^3/s]",       ## plot a histogram
     main="River Discharge",sub="log-normal model",
     xlim=c(0,2000),breaks=60,prob=TRUE)
lines(x,dlnorm(x,log.mu,log.sd),col="red") ## add the log-normal distribution with estimated parameters


### Extremal Types Theorem
##########################

### generate samples of random variables
N <- 100000                    ## length
x <- rnorm(N,mean=0,sd=1)      ## standard normal
## alterate the probability distribution by commenting out the former line
## and uncommenting one of the following
#x <- rlnorm(N,mean=0,sd=1)     ## standard log-normal
#x <- runif(N,min=0,max=1)      ## uniform
#x <- rexp(N,rate=1)            ## exponential

## plot sample
plot(x)

## plot histogram
hist(x,breaks=100,prob=TRUE)

## draw maxima with block size n
n <- 100
Mn <- apply(matrix(x,nrow=n),2,max)

## plot maxima M_n
plot(Mn)

## histogram of M_n with density estimate
hist(Mn,prob=TRUE)
lines(density(Mn),col=2)