Search code examples
pythonnumpyscipystatisticsdistribution

Calculating Return Value using Generalised Extreme Value Distribution (GEV)


I have annual temperature data of 30 years and want to calculate Return Value of this data using GEV distribution for 50 and 100 year Return Period.

My 30 year data:

data=[28.01,29.07,28.67,21.57,21.66,24.62,21.45,28.51,22.65,21.57,20.89,20.96,21.05,22.29,20.81,21.08,20.77,23.18,22.98,21.88,21.07,20.74,22.69,22.42,31.81,25.78,29.09,28.11,22.18,21.6]

How to find return value using GEV?


Solution

  • To estimate the return level of a given return period T, first estimate the parameters of the generalized extreme value distribution, and then compute the inverse of the survival function at 1/T of the fitted distribution. (The survival function SF(x) is just 1 - CDF(x). If you read about computing return levels, you'll typically see the problem stated as solving CDF(x) = 1 - 1/T. That is the same as solving SF(x) = 1/T.)

    Here's a script that uses scipy.stats.genextreme to estimate the return levels for your data at several return periods. The method genextreme.isf() is the inverse of the survival function.

    import numpy as np
    from scipy.stats import genextreme
    
    
    data = np.array([28.01, 29.07, 28.67, 21.57, 21.66, 24.62, 21.45, 28.51,
                     22.65, 21.57, 20.89, 20.96, 21.05, 22.29, 20.81, 21.08,
                     20.77, 23.18, 22.98, 21.88, 21.07, 20.74, 22.69, 22.42,
                     31.81, 25.78, 29.09, 28.11, 22.18, 21.6])
    
    # Fit the generalized extreme value distribution to the data.
    shape, loc, scale = genextreme.fit(data)
    print("Fit parameters:")
    print(f"  shape: {shape:.4f}")
    print(f"  loc:   {loc:.4f}")
    print(f"  scale: {scale:.4f}")
    print()
    
    # Compute the return levels for several return periods.
    return_periods = np.array([5, 10, 20, 50, 100])
    return_levels = genextreme.isf(1/return_periods, shape, loc, scale)
    
    print("Return levels:")
    print()
    print("Period    Level")
    print("(years)   (temp)")
    
    for period, level in zip(return_periods, return_levels):
        print(f'{period:4.0f}  {level:9.2f}')
    

    Output:

    Fit parameters:
      shape: -0.9609
      loc:   21.5205
      scale: 1.0533
    
    Return levels:
    
    Period    Level
    (years)   (temp)
       5      25.06
      10      29.95
      20      39.45
      50      67.00
     100     111.53