I created a synthetic series that is supposed to simulate a series of peak discharges in blocks of years in arid catchments. The magnitudes were simulated via the Lnorm dist.:
meanlog <- 1.86; sdlog = 1.63
x.lnorm <- rlnorm(1000000,meanlog = meanlog, sdlog = sdlog)
The occurances per year were simulated by the Negative Binomial dist.:
size = 4.72; mu = 2.2
x.binom <- rnbinom(100000,size=size,mu = mu)
Then I create the series like this:
flow <- list()
for (i in 1:length(x.binom)){
if (x.binom[i]!=0){
flow[[i]] <- sample(x.lnorm,x.binom[i])
}
else{
flow[[i]] <- 0
}
}
So I have a simulation of 100000 years of floods. I order to create the Block Maxima (AMS) series I extract the largest flood from every year in the series. Next, I fit the GEV to the series with both the methods: MLE and PWM.
Now, I get this result for the MLE: and this result for the PWM:
Both results are not perfect but of course the MLE method doesn't fit at all. My first question is: Why does the MLE estimation render so far from the simulated series at ARI (average recurrence intervals) that are above 10 and the pwm fit is relatively good?
I have run experiments of using larger params in the N.Binom dist. so the block maxima have much more events per block (complies better with the GEV theory) and the fits were (as expected) much better for both estimation methods. But this only left me even more puzzled as to why the pwm does manage to emulate the data and the mle does not. What is causing the shape parameter to be so large with the mle method?
Just for clarification: all the zeros are removed from the AMS (so this is not the reason for the mle shooting up) and the probability calc. is:
# Zero flow years treatment
fracZero <- sum(ams!=0)/length(ams)
amsnozero <- ams[ams!=0]
Tvec <- c(seq(from = 10, to = 1.1, by = -0.1) %o% 10^(2:0))
Pvec <- 1-(1./Tvec) #nonExceedance prob
Pvec <- (Pvec-1+fracZero)/fracZero; Pvec[Pvec<0] <- 0 #from Haan 2002