### Name: GevModelling
### Title: Generalized Extreme Value Modelling
### Aliases: dgev pgev qgev rgev devd pevd qevd revd GevFit gevSim gevFit
###   print.gevFit plot.gevFit summary.gevFit gevrlevelPlot hillPlot
###   shaparmPlot shaparmPickands shaparmHill shaparmDEHaan
### Keywords: models

### ** Examples

## SOURCE("fExtremes.52A-GevModelling")

## *gev  -
   xmpExtremes("\nStart: GEV Frechet >")
   # Create and plot 1000 GEV/Frechet distributed rdv:
   par(mfrow = c(3, 3))
   r = rgev(n = 1000, xi = 1)
   plot(r, type = "l", main = "GEV/Frechet Series")
   ## Plot empirical density and compare with true density:
   ## Omit values greater than 500 from plot
   hist(r[r<10], n = 25, probability = TRUE, xlab = "r", 
     xlim = c(-5, 5), ylim = c(0, 1.1), main = "Density")
   x = seq(-5, 5, by=0.01)
   lines(x, dgev(x, xi = 1), col = 2)
   ## Plot df and compare with true df:
   plot(sort(r), (1:length(r)/length(r)), 
     xlim = c(-3, 6), ylim = c(0, 1.1),
     cex = 0.5, ylab = "p", xlab = "q", main = "Probability")
   q = seq(-5,5, by=0.1)
   lines(q, pgev(q, xi=1), col=2)
## Compute quantiles, a test:
   qgev(pgev(seq(-5, 5, 0.25), xi = 1), xi = 1) 
   
## *gev -
   xmpExtremes("\nNext: GEV Gumbel >")
   # Create and plot 1000 Gumbel distributed rdv:
   ##> r = rgev(n = 1000, xi = 0)
   ##> plot(r, type = "l", main = "Gumbel Series")
   ## Plot empirical density and compare with true density:
   ##>hist(r[abs(r)<10], nclass = 25, freq = FALSE, xlab = "r", 
   ##>   xlim = c(-5,5), ylim = c(0,1.1), main = "Density")
   ##>x = seq(-5, 5, by = 0.01)
   ##>lines(x, dgev(x, xi = 0), col=2)
   ## Plot df and compare with true df:
   ##>plot(sort(r), (1:length(r)/length(r)), 
   ##>  xlim = c(-3, 6), ylim = c(0, 1.1),
   ##>   cex=0.5, ylab = "p", xlab="q", main="Probability")
   ##>q = seq(-5, 5, by = 0.1)
   ##>lines(q, pgev(q, xi = 0), col = 2)
   ## Compute quantiles, a test:
   ##>qgev(pgev(seq(-5, 5, 0.25), xi = 0), xi = 0)   

## *gev -
   xmpExtremes("\nNext: GEV Weibull >")
   # Create and plot 1000 Weibull distributed rdv:
   r = rgev(n = 1000, xi = -1)
   plot(r, type = "l", main = "Weibull Series")
   ## Plot empirical density and compare with true density:
   hist(r[abs(r)<10], nclass = 25, freq = FALSE, xlab = "r", 
     xlim=c(-5,5), ylim=c(0,1.1), main="Density")
   x = seq(-5, 5, by=0.01)
   lines(x, dgev(x, xi = -1), col = 2)
   ## Plot df and compare with true df:
   plot(sort(r), (1:length(r)/length(r)), 
     xlim = c(-3, 6), ylim = c(0, 1.1),
     cex = 0.5, ylab = "p", xlab = "q", main = "Probability")
   q=seq(-5, 5, by = 0.1)
   lines(q, pgev(q, xi = -1), col = 2)
   ## Compute quantiles, a test:
   qgev(pgev(seq(-5, 5, 0.25), xi = -1), xi = -1)   

## gevSim -
## gevFit -
   # Simulate GEV Data:
   xmpExtremes("\nStart: Simulte GEV Sample >")
   # Use default length n=1000
   x = gevSim(model = list(shape = 0.25, location =0 , scale = 1))
   # Fit GEV Data by Probability Weighted Moments:
   fit = gevFit(x, type = "pwm") 
   print(fit)
   # Summarize Results:
   par(mfcol = c(3, 2))
   summary(fit)
   
## gevFit -
   # Fit GEV Data by Max Log Likelihood Method:
   xmpExtremes("\nNext: Estimate Parameters >")
   fit = gevFit(x, type = "mle") 
   print(fit)
   # Summarize Results:
   summary(fit) 
    
## gevSim -
## gevFit -
   # Simulate Gumbel Data:
   xmpExtremes("\nNext: Simulte Gumbel Sample >")
   # Use default length n=1000
   ##> x = gevSim(model = list(shape = 0, location = 0, scale = 1))
   # Fit Gumbel Data by Probability Weighted Moments:
   ##> fit = gevFit(x, type = "pwm", gumbel = TRUE) 
   ##> print(fit)
   # Summarize Results:
   ##> par(mfcol = c(3, 2))
   ##> summary(fit)  
   
## Fit Gumbel Data by Max Log Likelihood Method:
   xmpExtremes("\nNext: Estimate Parameters >")
   ##> fit = gevFit(x, type = "mle", gumbel = TRUE) 
   ##> print(fit)
   # Summarize Results:
   ##> summary(fit)     
   ##> xmpExtremes("Press any key to continue >") 
   
## Return levels based on GEV Fit:
   # BMW Stock Data:
   xmpExtremes("\nNext: Compute BMW Return Levels >")
   par(mfrow = c(2, 1))
   data(bmw)
   # Fit GEV to monthly Block Maxima:
   fit = gevFit(-bmw, block = "month") 
   # Calculate the 40 month return level
   gevrlevelPlot(fit, k.block = 40, main = "BMW: Return Levels")   
   
## Return levels based on GEV Fit:
   xmpExtremes("\nNext: Compute SIEMENS Return Levels >")
   # Siemens Stock Data:
   data(siemens)
   # Fit GEV to monthly Block Maxima:
   fit = gevFit(-siemens, block = "month") 
   # Calculate the 40 month return level
   gevrlevelPlot(fit, k.block = 40, main = "SIEMENS: Return Levels")
   
## Interactive Plot:
   ##> par(mfrow = c(1, 1), ask = TRUE)
   ##> plot(fit)
   
## hillPlot -
   xmpExtremes("\nStart: Hill Estimator >")
   # Hill plot of heavy-tailed Danish fire insurance data 
   # and BMW stock data for estimated 0.999 quantile
   par(mfrow = c(2, 2))
   data(bmw)
   hillPlot(bmw)
   hillPlot(bmw, option = "quantile", end = 500, p = 0.999)
   data(danish)
   hillPlot(danish)
   hillPlot(danish, option = "quantile", end = 500, p = 0.999)
   
## shaparmPlot -
   xmpExtremes("\nNext: Shape Parameter Plots >")
   par(mfcol = c(3, 2), cex = 0.6)
   data(bmw)
   shaparmPlot(bmw)

## shaparmPlot -
   xmpExtremes("\nNext: Simulated Frechet Data >")
   par(mfcol = c(3, 2), cex = 0.6)
   set.seed(4711)
   x = rgev(10000, xi = 1/4)
   shaparmPlot(x, revert = TRUE, both.tails = FALSE)
   lines(c(0.01, 0.1), c(4, 4), col = "steelblue3") # True Value



