Performance of the generalized least-squares method for the Gumbel distribution and its application to annual maximum wind speeds

J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132

Contents lists available at SciVerse ScienceDirect

Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia

Performance of the generalized least-squares method for the Gumbel distribution and its application to annual maximum wind speeds H.P. Hong n, S.H. Li, T.G. Mara Department of Civil and Environmental Engineering, University of Western Ontario, London, ON, Canada N6A 5B9

art ic l e i nf o

a b s t r a c t

Article history: Received 16 July 2012 Received in revised form 8 May 2013 Accepted 19 May 2013 Available online 28 June 2013

The Gumbel model is often used to fit annual maximum wind speed or wind velocity pressure. The commonly used fitting methods include the method of moments, the method of maximum likelihood, the method of L-moments, and the Lieblein BLUE (i.e., generalized least-squares method (GLSM)). Previously, the coefficients of the estimators for the latter method have not been available for large sample size, and the relative performance of the GLSM to other fitting methods such as the method of L-moments is unknown. In this study, we evaluate these coefficients for a sample size up to 100, and identify trends in the calculated coefficients. The relative performance of commonly used fitting methods for the Gumbel distribution, including the GLSM, is evaluated in terms of efficiency, bias, and root-mean square error. We illustrate their application and impact on the estimated return period values of the annual maximum wind speed for 14 locations in Canada. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Distribution fitting method Generalized least-squares method Lieblein BLUE Gumbel Annual maximum wind speed

1. Introduction In extreme value theory, the linear-stability postulate leads to only three types of extreme value probability distributions (Leadbetter et al., 1983; Castillo, 1988). One of the three types, commonly referred to as the Gumbel distribution (or extreme value type I), is often used to model the annual maxima of environmental phenomena. For example, this distribution has been used to model the annual maximum wind speed by Yip et al. (1995) for Canada, by Peterka and Shahid (1998) for the United States, by Frank (2001) for Denmark, and by Sacré (2002) for France. The Gumbel model has also been considered in the literature for the annual maximum wind velocity pressure, which is directly proportional to squared wind speed. This could be justified based on the extremal distribution of an independent and identically distributed Weibull sequence (Hong, 1994; Cook and Harris, 2004). The performance of fitting methods such as the method of moments (MOM), the method of maximum likelihood (MML), and the probability weighted moments (PWM) for the Gumbel model is compared by Phien (1987) in terms of the root mean square error (RMSE) and efficiency, which is defined based on the ratio of the variance of the predicted quantile of the models. It was concluded that the PWM method is best if one is interested in unbiased estimates, and that the MML method is the most

n

Corresponding author. Tel.: +1 519 661 2111x88315; fax: +1 519 661 3779. E-mail address: [email protected] (H.P. Hong).

0167-6105/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jweia.2013.05.012

appropriate based on the minimum RMSE criterion. Both the MML and PWM outperform the MOM. Moreover, Hosking (1990) showed that the estimates obtained from the method of the L-moments (MLM) are more accurate than the maximum likelihood estimates if the sample size is small, and that the procedures based on the PWM and on the MLM are equivalent. The term PWM will be dropped and replaced by the method of L-moments (MLM) for the remaining discussion of this study. In addition to the methods mentioned above, the Lieblein BLUE (best linear unbiased estimators) (Lieblein 1974) is also used for Gumbel distribution fitting. The method is an application of the generalized least-squares method (GLSM) (Draper and Smith, 1998) to the Gumbel model; the use of the GLSM was considered by Lloyd (1952) to estimate the scale and location parameters of probability distributions. Harris (1996) indicated that the coefficients of the Lieblein BLUE (i.e., GLSM) for sample size greater than 24 are not available in the literature, and proposed an approximation to the Lieblein BLUE by retaining the variance and neglecting the covariance of the order statistics in the regression analysis. Harris indicated that this approximation is a viable alternative to the complete Lieblein BLUE for dealing with annual maxima. Moreover, he suggested that the probability distribution of the standardized 50-year return period value is slightly skewed and insensitive to the number of years of data employed for the distribution fitting. However, a comparison among the GLSM, the MML and the MLM is lacking in the literature. The main objectives of the present study are to compare the performance of the MLM and GLSM for the Gumbel model, and to investigate the impact associated with the choice of a fitting

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132\94\maccounttest15=tx1

Table 1 Coefficients of the BLUE (for n¼ 5 to 15, 20, 30, 40, 50, 100). cu,i n ¼5 0.418934 0.246282 0.167609 0.108824 0.058350 n ¼6 0.355450 0.225488 0.165620 0.121054 0.083522 0.048867 n ¼7 0.309007 0.206260 0.158590 0.123223 0.093747 0.067331 0.041841 n ¼8 0.273535 0.189428 0.150200 0.121175 0.097142 0.075904 0.056132 0.036485 n ¼9 0.245539 0.174882 0.141789 0.117357 0.097218 0.079569 0.063400 0.047957 0.032291 n ¼10 0.222867 0.162308 0.133845 0.112868 0.095636 0.080618 0.066988 0.054193 0.020485 0.019660 0.018848 0.018049 0.017259 0.016478 0.015703 0.014933 0.014165 0.013398 0.012630 0.011857 0.011076 0.010283 0.009472 0.008633 0.007748 0.006771 n ¼50 0.049601 0.043109 0.039894 0.037559 0.035680 0.034085 0.032688

ca,i

−0.503127 0.006534 0.130455 0.181656 0.184483 −0.459273 −0.035992 0.073199 0.126724 0.149534 0.145807 −0.423700 −0.060698 0.036192 0.087339 0.114868 0.125859 0.120141 −0.394187 −0.075767 0.011125 0.058928 0.087162 0.102728 0.108074 0.101937 −0.369242 −0.085203 −0.006486 0.037977 0.065574 0.082655 0.091965 0.094369 0.088391 −0.347830 −0.091158 −0.019210 0.022179 0.048671 0.066065 0.077021 0.082771 0.015501 0.016186 0.016790 0.017318 0.017772 0.018155 0.018470 0.018717 0.018897 0.019009 0.019052 0.019021 0.018912 0.018717 0.018422 0.018003 0.017417 0.016539 −0.124201 −0.066835 −0.049770 −0.039020 −0.031279 −0.025300 −0.020478

cu,i 0.041748 0.028929 n¼ 11 0.204123 0.151385 0.126522 0.108226 0.093234 0.080222 0.068485 0.057578 0.047159 0.036886 0.026180 n¼ 12 0.188361 0.141833 0.119838 0.103673 0.090455 0.079018 0.068747 0.059266 0.050303 0.041628 0.032984 0.023894 n¼ 13 0.174916 0.133422 0.113759 0.099323 0.087540 0.077368 0.068265 0.059900 0.052048 0.044528 0.037177 0.029790 0.021965 n¼ 14 0.163309 0.125966 0.108230 0.095223 0.084619 0.075484 0.067331 0.018410 0.017850 0.017300 0.016759 0.016227 0.015703 0.015185 0.014674 0.014167 0.013665 0.013167 0.012672 0.012180 0.011688 0.011198 0.010708 0.010216 0.009723 0.009225 0.008723 0.008213 0.007694 0.007161 0.006607 0.006022 0.005374

ca,i 0.083551 0.077940 −0.329210 −0.094869 −0.028604 0.010032 0.035284 0.052464 0.064071 0.071381 0.074977 0.074830 0.069644 −0.312840 −0.097086 −0.035655 0.000534 0.024548 0.041278 0.053053 0.061113 0.066122 0.068357 0.067670 0.062906 −0.298313 −0.098284 −0.041013 −0.006997 0.015836 0.032014 0.043710 0.052101 0.057862 0.061355 0.062699 0.061699 0.057329 −0.285316 −0.098775 −0.045120 −0.013039 0.008690 0.024282 0.035768 0.010339 0.010952 0.011513 0.012025 0.012493 0.012917 0.013300 0.013644 0.013950 0.014219 0.014454 0.014653 0.014818 0.014949 0.015046 0.015108 0.015134 0.015124 0.015075 0.014984 0.014848 0.014661 0.014413 0.014089 0.013661 0.013048

cu,i 0.059866 0.052891 0.046260 0.039847 0.033526 0.027131 0.020317 n¼ 15 0.153184 0.119314 0.103196 0.091384 0.081767 0.073495 0.066128 0.059401 0.053140 0.047217 0.041529 0.035984 0.030482 0.024887 0.018894 n¼ 20 0.117208 0.094565 0.083697 0.075756 0.069317 0.063808 0.058934 0.054520 0.050452 0.046652 0.043062 0.039637 0.036341 0.033140 0.030013 0.026920 0.023829 0.020710 0.017477 0.013961 n¼ 30 0.080216 0.067253 0.060943 0.056372 0.015237 0.014969 0.014711 0.014463 0.014222 0.013990 0.013765 0.013546 0.013333 0.013126 0.012923 0.012726 0.012533 0.012345 0.012160 0.011979 0.011801 0.011627 0.011456 0.011288 0.011122 0.010959 0.010799 0.010640 0.010485 0.010331

ca,i 0.044262 0.050418 0.054624 0.057083 0.057829 0.056652 0.052642 −0.273606 −0.098767 −0.048285 −0.017934 0.002773 0.017779 0.028988 0.037452 0.043798 0.048415 0.051534 0.053267 0.053604 0.052334 0.048648 −0.228751 −0.095077 −0.055998 −0.032080 −0.015396 −0.002961 0.006674 0.014319 0.020471 0.025448 0.029473 0.032695 0.035221 0.037120 0.038447 0.039204 0.039380 0.038936 0.037693 0.035183 −0.175813 −0.083973 −0.056827 −0.039960 −0.006825 −0.006046 −0.005321 −0.004644 −0.004010 −0.003416 −0.002857 −0.002330 −0.001832 −0.001361 −0.000915 −0.000493 −0.000091 0.000291 0.000655 0.001001 0.001332 0.001647 0.001948 0.002236 0.002511 0.002774 0.003026 0.003267 0.003498 0.003719

cu,i

ca,i

0.052645 0.049496 0.046719 0.044219 0.041934 0.039817 0.037836 0.035967 0.034192 0.032495 0.030865 0.029292 0.027767 0.026283 0.024833 0.023411 0.022013 0.020630 0.019260 0.017892 0.016526 0.015147 0.013746 0.012306 0.010795 0.009132 n¼ 40 0.061215 0.052443 0.048135 0.045004 0.042480 0.040336 0.038454 0.036765 0.035226 0.033804 0.032479 0.031234 0.030057 0.028937 0.027866 0.026838 0.025849 0.024892 0.023965 0.023063 0.022184 0.021326 0.007827 0.007698 0.007570 0.007444 0.007316 0.007190 0.007065 0.006940 0.006816 0.006692 0.006568 0.006445 0.006322 0.006199 0.006077 0.005954 0.005832 0.005710 0.005587 0.005465 0.005343 0.005220 0.005097 0.004974 0.004850 0.004726

−0.027950 −0.018822 −0.011566 −0.005638 −0.000691 0.003491 0.007066 0.010146 0.012811 0.015123 0.017132 0.018872 0.020375 0.021662 0.022753 0.023660 0.024395 0.024977 0.025352 0.025624 0.025694 0.025594 0.025296 0.024763 0.023917 0.022539 −0.144857 −0.074407 −0.053497 −0.040387 −0.030993 −0.023775 −0.017984 −0.013200 −0.009162 −0.005700 −0.002695 −0.000062 0.002263 0.004327 0.006169 0.007818 0.009298 0.010627 0.011823 0.012899 0.013864 0.014729 0.006421 0.006515 0.006604 0.006690 0.006768 0.006843 0.006914 0.006980 0.007041 0.007098 0.007151 0.007200 0.007245 0.007285 0.007321 0.007353 0.007381 0.007404 0.007423 0.007439 0.007449 0.007456 0.007457 0.007455 0.007447 0.007435

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132

123

Table 1 (continued ) cu,i

ca,i

cu,i

ca,i

cu,i

ca,i

cu,i

ca,i

0.031436 0.030296 0.029246 0.028269 0.027353 0.026489 0.025668 0.024886 0.024138 0.023418 0.022725 0.022055 0.021407 0.020776 0.020163 0.019565 0.018981

−0.016470 −0.013067 −0.010130 −0.007562 −0.005295 −0.003276 −0.001466 0.000165 0.001642 0.002984 0.004208 0.005326 0.006351 0.007290 0.008152 0.008944 0.009671

n¼ 100 0.025703 0.023143 0.021833 0.020884 0.020122 0.019477 0.018913 0.018409 0.017952 0.017533 0.017144 0.016780 0.016439 0.016115 0.015809 0.015517

−0.075640 −0.045356 −0.036322 −0.030574 −0.026392 −0.023129 −0.020470 −0.018237 −0.016320 −0.014646 −0.013166 −0.011843 −0.010650 −0.009568 −0.008575 −0.007666

0.010179 0.010029 0.009881 0.009735 0.009590 0.009447 0.009305 0.009165 0.009026 0.008889 0.008752 0.008617 0.008483 0.008350 0.008218 0.008086 0.007956

0.003931 0.004134 0.004328 0.004515 0.004693 0.004864 0.005028 0.005185 0.005336 0.005480 0.005617 0.005749 0.005875 0.005995 0.006109 0.006219 0.006323

0.004601 0.004475 0.004349 0.004221 0.004092 0.003962 0.003831 0.003697 0.003561 0.003422 0.003280 0.003133 0.002980 0.002818 0.002639

0.007418 0.007395 0.007368 0.007334 0.007295 0.007248 0.007195 0.007134 0.007065 0.006985 0.006895 0.006790 0.006667 0.006520 0.006327

method on the estimated return period values of the annual maximum wind speed for a number of sites in Canada. To achieve the first objective, the equations required to calculate the coefficients of the best linear unbiased estimators based on the GLSM method are programmed, and the coefficients for a sample size up to 100 are calculated. These coefficients are used to calibrate empirical equations which can be used as a crude approximation to calculate the coefficients of the BLUE for sample size greater than about 15. The MOM and MML are also included in the comparison for completeness. The impact of using these methods to estimate the return period values of the maximum wind speed is illustrated using wind records from 14 Canadian meteorological stations.

2. Model and methods 2.1. Gumbel model and the fitting methods The most commonly used probabilistic model to deal with the extreme wind speeds is the Gumbel model, which is also known as extreme value type I distribution. The Gumbel model is given by Castillo (1988), FðxÞ ¼ expð−expð−ðx−uÞ=aÞÞ

ð1Þ

where F(x) denotes the cumulative distribution function, x denotes the value of the random variable X, and u and a are the location and scale parameters, respectively. Y ¼(X−u)/a represents the standardized or reduced extreme variate. The mean of X,pμffiffiffiX, equals u þ aγ and the standard deviation of X, sX, equals aπ= 6, where γ≈0.5772 is the Euler constant. The quantile (or the T-year return period value) of X, xT, can be estimated using, xT ¼ u þ ayT

ð2Þ

where yT ¼ −lnð−lnð1−1=TÞÞ. If the MOM is considered, the model parameters are calculated using the sample mean m and sample standard deviation s replacing those of the population. If the MML is considered, the ^ are to be obtained by estimates of a and u, denoted by a^ and u, maximizing the log-likelihood function, L, (Phien, 1987). If the MLM is adopted for the distribution fitting, a^ ¼ ð2b1 −b0 Þ=ln 2 and u^ ¼ b0 −γ a^ , where bi is given by Hosking et al. (1985) and Hosking (1990), ðj−1Þðj−2Þ⋯ðj−iÞ xj:n ; for i ¼ 0; 1; 2; … j ¼ 1 ðn−1Þðn−2Þ⋯ðn−iÞ n

bi ¼ n−1 ∑

ð3Þ

in which xj:n denotes the the jth ordered sample (in ascending order) of a set of random samples of size n. For the Gumbel model, the order statistics xi:n can be expressed as, xi:n ¼ u þ ayi:n

ð4Þ

where yi:n ¼ −lnð−lnðFðxi:n ÞÞÞ represents the reduced order statistics. Equations for the probability density function of yi:n, the moments of yi:n, and the second order moments of yi:n and yj:n are given in Lieblein (1953, 1974). Given the samples, the model parameters u and a can be estimated using the ordinary least-squares method if the variance–covariance matrix of yi:n and yj:n equals the identity matrix except for a scaling constant. As shown by Lieblein (1953, 1974) the variance–covariance matrix of the order statistics for the considered model does not satisfy such a criterion, and the ordinary least-squares method should not be applied. Consequently, the GLSM should be considered (Lloyd, 1952; Draper and Smith, 1998). In such a case, the minimization of the sum of the square of residuals (or difference) between the left hand side and the right hand side of Eq. (4), Ve, given by Lloyd (1952), V e ¼ ðx−AθÞT B−1 ðx−AθÞ

ð5Þ

results in, θ ¼ ðc u

c a ÞT x

ð6Þ

where ðc u

ca ÞT ¼ ðAT B−1 AÞ−1 AT B−1 ;

ð7Þ

cu and ca are both 1  n matrices (with elements denoted by cu,i and ca,i, respectively),xT ¼ ½x1:n ; …; xn:n , IT ¼ ½1; …; 1, θT ¼ ½u; a, A ¼ ½I α, αT ¼ ½α1:n ; …; αn:n , and B ¼ Covarðyi:n ; yj:n Þ represents the variance–covariance matrix. The equations to evaluate the elements in α and B (i.e., αki:n ¼ Eðyki:n Þ and βi;j:n ¼ EðY i:n Y j:n Þ−αi:n αj:n ) were given by Lieblein (1953, 1974) (see Appendix). The Lieblein approach was implemented, and it was concluded that it is only suitable to evaluate α and B for n up to 20 (due to alternate signs in the formulation). Solution of Eq. (7) provides the weight, or coefficients, for using the ordered samples to estimate the location and scale parameters of the Gumbel model for the Lieblein method (see Eq. (6)), that can be re-written as, n

n

i¼1

i¼1

u ¼ ∑ cu;i xi:n ; and a ¼ ∑ ca;i xi:n

ð8Þ

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132\94\maccounttest15=tx1 Table 2 Approximate BLUE coefficients estimated by using Eqs. (11a) and (11b) for n equal to 15, 30 and 50. cu,i

Fig. 1. Variation of the ratio ru(z).

n¼ 15 0.162 0.119 0.102 0.090 0.081 0.073 0.066 0.059 0.053 0.047 0.041 0.036 0.030 0.024 0.018 n¼ 30 0.082 0.067 0.061 0.056 0.053 0.049 0.047 0.044

ca,i

−0.277 −0.075 −0.036 −0.015 0.000 0.010 0.018 0.024 0.029 0.032 0.034 0.034 0.034 0.033 0.031 −0.173 −0.072 −0.048 −0.034 −0.024 −0.016 −0.010 −0.005

cu,i

ca,i

cu,i

ca,i

cu,i

ca,i

0.042 0.040 0.038 0.036 0.034 0.032 0.031 0.029 0.028 0.026 0.025 0.023 0.022 0.021 0.019 0.018 0.016 0.015 0.014 0.012 0.011 0.009 n¼ 50 0.050 0.043

−0.001 0.002 0.005 0.008 0.010 0.013 0.014 0.016 0.017 0.018 0.019 0.020 0.021 0.021 0.021 0.021 0.021 0.021 0.021 0.021 0.020 0.019

0.040 0.038 0.036 0.034 0.033 0.031 0.030 0.029 0.028 0.027 0.026 0.026 0.025 0.024 0.023 0.023 0.022 0.021 0.021 0.020 0.020 0.019 0.018 0.018 0.017

−0.047 −0.037 −0.030 −0.025 −0.020 −0.016 −0.013 −0.010 −0.008 −0.006 −0.004 −0.002 0.000 0.001 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.010 0.011 0.012

0.017 0.016 0.016 0.015 0.015 0.014 0.014 0.013 0.013 0.012 0.012 0.011 0.011 0.010 0.010 0.009 0.009 0.008 0.008 0.007 0.007 0.006 0.006

0.012 0.013 0.013 0.014 0.014 0.014 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.014 0.014 0.014

−0.123 −0.064

Fig. 2. Variation of the ratio ra(z).

Fig. 3. Coefficients for ca,1.

The coefficients cu,i and ca,i, which are known as the coefficients of the best linear unbiased estimators (BLUE) (i.e., the estimators given by Eq. (7) having minimum variance), are tabulated in Lieblein (1974) for n≤16 and to 6 decimal points, and in Cook (1985) for n≤24 to 3 decimal points. For n≤30, one could infer the coefficients from Balakrishnan and Chan (1992b). For n greater than 30, coefficients of the BLUE are not available in the literature. The numerical difficulty in evaluating αi:n and βi,j:n for large n values is discussed by Balakrishnan and Chan (1992a) and by Harris (1996).

2.2. Estimated coefficients for the BLUE Estimation of cu and ca for n up to 100 is carried out by using the implemented procedure discussed in the Appendix, although our computing implementation can be used to calculate the coefficients for larger values of n. The coefficients are tabulated in Table 1 for n equal to 5 to 15, 20, 30, 40, 50, and 100. The values are in agreement with those reported by Lieblein (1974) for up to 6 decimal points for n≤16 (except in a few values for up to 5 decimal points), and show good agreement with those reported by Cook (1985) for n≤24. The conditions that the calculated cu,i and

Fig. 4. Comparison of the coefficients by including and excluding the covariance: (a) comparison of cu,i for a few selected n values; (b) comparison of ca,i for a few selected n values.

ca,i satisfy ∑ni¼ 1 cu;i ¼ 1 and ∑ni¼ 1 ca;i ¼ 0 (Lieblein, 1953; Harris, 1996) are verified. To appreciate the behaviour of the coefficients, and to potentially develop empirical equations to estimate the coefficients of the BLUE, the ratio ru(z) (ru(z) ¼cu,i/cu,1) is plotted in Fig. 1 while the ratio ra(z) (ra(z) ¼ca,i/ca,1) is depicted in Fig. 2. The abscissa in

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132

each figure represents z ¼(i−1)/(n−1). Figs. 1 and 2 show that ra(z) equals zero at about 0.288 and that rc(z) and ra(z) are well behaved. Through trial and error, it was found that each of the ratios rc(z) and ra(z) could be crudely approximated by a single curve for n≥15. In fact, for n≥15, it was found that the following simple empirical equations, r u ðzÞ ¼ expð−1:386z0:571 Þ−0:137z2:269 ;

ð9aÞ

and, r a ðzÞ ¼ expð−4:66z0:507 Þ−0:907atanðzÞ þ 0:591z;

ð9bÞ

fit the curves well, where z¼ (i−1)/(n−1). The fitting is carried out by minimizing sum of the squared errors. Moreover, the values of ca,1 which are shown in Fig. 3, can be approximated by ca(n), ca ðnÞ ¼ 0:2542lnðnÞ−2:22  10−5 n2 −0:0058n−0:8832

ð10Þ

The absolute value of the error (i.e., the difference between the exact value and the predicted value by the approximations) by using Eqs. (9a), (9b) and (10) is less than 0.04, 0.09 and 0.01, respectively, for n ranging from 15 to 100. With these equations, and considering that ∑ni¼ 1 cu;i ¼ 1, the coefficients of the BLUE may be approximated by, cu;i ¼ r u ðzÞ=C s

ð11aÞ

and ca;i ¼ ca ðnÞ  r a ðzÞ; where

ð11bÞ

125

The approximate BLUE coefficients estimated by Eqs. (11a) and (11b) for selected values of n are shown in Table 2. Comparison of the results shown in Tables 1 and 2 indicate that the approximate coefficients are close to the “exact” values presented in Table 1. It is emphasized that the provided approximations are to be used for crude and preliminary analysis only. Harris (1996) advocated the use of the inverse of the variance of order statistics as weights and neglecting the covariance for the calculation of sets of coefficients for the distribution fitting due to the unavailability of coefficients of the BLUE. This approximation, explained in terms of GLSM, is replacing the variance– covariance matrix B shown in Eqs. (5)–(7) by the corresponding variance matrix. This approximation will be herein referred to as GLSM-V. The calculated coefficients based on variance alone are shown in Fig. 4 and are compared with those of the BLUE, indicating that the former and the latter do not follow a common trend. However, it will be seen in the next section that the distribution fitting obtained by using the GLSM-V is comparable to that of the GLSM.

3. Comparison of the performance The performance of the fitting methods for the Gumbel model can be evaluated in terms of the relative efficiency, Bias and RMSE (Phien, 1987). The relative efficiency of method I and method II in estimating a variable ζ, representing u or a, or xT, is defined by, ^ ^ Ef f ¼ varianceðζjmethodIIÞ=varianceð ζjmethodIÞ;

C s ¼ ∑nj¼ 1 r u ðzÞ.

ð12Þ

Table 3 Variance, Bias and RMSE for different methods (Bias and RMSE are defined in terms of normalized values as shown in Eqs. (13) and (14), except for u which are calculated based on relative values). Parameter

n ¼10 a u X30 X50 x100 x500 n ¼20 a u X30 X50 x100 x500 n ¼30 a u X30 X50 x100 x500 n ¼50 a u X30 X50 x100 x500 n ¼100 a u X30 X50 x100 x500

Variance (  10)

Bias (  100)

RMSE (  10)

MOM

MML

MLM

GLSM

GLSM-V

MOM

MML

MLM

GLSM

GLSM-V

MOM

MML

MLM

GLSM

GLSM-V

0.89 1.17 12.43 15.95 21.46 37.50

0.60 1.15 9.65 12.17 16.08 27.37

0.86 1.14 12.26 15.69 21.04 36.61

0.71 1.13 10.81 13.74 18.29 31.46

0.74 1.15 11.10 14.13 18.84 32.49

−5.06 2.58 −4.30 −4.40 −4.50 −4.64

−7.95 3.72 −6.86 −7.00 −7.15 −7.36

−0.32 −0.15 −0.37 −0.36 −0.36 −0.35

−0.24 −0.17 −0.29 −0.29 −0.28 −0.27

−0.22 −0.15 −0.26 −0.25 −0.25 −0.24

3.03 3.43 3.32 3.27 3.22 3.15

2.58 3.41 2.98 2.91 2.85 2.76

2.93 3.37 3.27 3.21 3.15 3.08

2.67 3.37 3.07 3.00 2.94 2.85

2.73 3.39 3.11 3.05 2.98 2.90

0.49 0.58 6.66 8.58 11.59 20.36

0.31 0.57 4.90 6.18 8.17 13.91

0.42 0.56 6.09 7.78 10.42 18.09

0.33 0.56 5.19 6.57 8.71 14.91

0.36 0.57 5.45 6.92 9.21 15.84

−2.46 1.33 −2.07 −2.12 −2.17 −2.25

−3.84 1.79 −3.31 −3.38 −3.45 −3.55

0.05 −0.12 0.01 0.02 0.02 0.03

0.02 −0.12 −0.01 −0.01 0.00 0.00

0.04 −0.13 0.00 0.01 0.01 0.02

2.23 2.42 2.42 2.38 2.35 2.31

1.79 2.39 2.10 2.04 2.00 1.93

2.05 2.38 2.31 2.26 2.22 2.16

1.82 2.37 2.13 2.08 2.03 1.96

1.89 2.40 2.18 2.13 2.09 2.03

0.34 0.39 4.51 5.82 7.87 13.87

0.20 0.37 3.25 4.10 5.41 9.21

0.27 0.37 4.00 5.11 6.84 11.86

0.21 0.37 3.38 4.27 5.66 9.66

0.24 0.38 3.61 4.58 6.09 10.47

−1.66 0.87 −1.41 −1.44 −1.48 −1.52

−2.61 1.17 −2.26 −2.31 −2.35 −2.42

0.01 −0.10 −0.02 −0.01 −0.01 0.00

−0.04 −0.09 −0.07 −0.06 −0.06 −0.06

−0.02 −0.11 −0.05 −0.04 −0.04 −0.03

1.84 1.97 1.99 1.96 1.93 1.90

1.45 1.93 1.70 1.66 1.62 1.56

1.66 1.93 1.87 1.83 1.80 1.75

1.46 1.92 1.72 1.67 1.64 1.58

1.54 1.95 1.77 1.73 1.70 1.65

0.21 0.23 2.77 3.58 4.85 8.55

0.12 0.22 1.94 2.45 3.24 5.51

0.16 0.22 2.39 3.06 4.09 7.09

0.13 0.22 1.99 2.51 3.33 5.67

0.14 0.23 2.16 2.75 3.65 6.28

−1.02 0.54 −0.86 −0.88 −0.90 −0.93

−1.52 0.68 −1.32 −1.35 −1.38 −1.41

0.04 −0.07 0.02 0.02 0.03 0.03

0.02 −0.07 0.00 0.00 0.01 0.01

0.04 −0.08 0.02 0.02 0.02 0.03

1.45 1.52 1.56 1.54 1.52 1.49

1.11 1.49 1.31 1.28 1.24 1.20

1.28 1.49 1.45 1.42 1.39 1.35

1.12 1.49 1.32 1.28 1.25 1.21

1.19 1.51 1.37 1.34 1.31 1.28

0.11 0.12 1.41 1.82 2.47 4.37

0.06 0.11 0.97 1.22 1.62 2.75

0.08 0.11 1.19 1.51 2.02 3.50

0.06 0.11 0.98 1.24 1.64 2.79

0.07 0.11 1.07 1.36 1.81 3.12

−0.52 0.37 −0.41 −0.42 −0.44 −0.46

−0.73 0.43 −0.60 −0.62 −0.64 −0.66

0.01 0.06 0.03 0.03 0.03 0.02

0.04 0.06 0.06 0.06 0.06 0.05

0.02 0.07 0.04 0.04 0.04 0.03

1.03 1.08 1.11 1.09 1.08 1.06

0.78 1.06 0.92 0.90 0.88 0.85

0.90 1.06 1.02 1.00 0.98 0.95

0.78 1.05 0.93 0.90 0.88 0.85

0.84 1.07 0.97 0.95 0.93 0.90

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132\94\maccounttest15=tx1

Table 4 Selected Canadian meteorological stations and estimated return period values of hourly-mean wind speed (km/h) by considering that the annual maximum wind speed is a Gumbel variate. Location

Province

Victoria Int’l A Whitehorse A Yellowknife A Iqaluit A Edmonton Int’l A Regina Int’l A Winnipeg Int’l A Ottawa Int’l A Toronto Int’l A Quebec Int’l A Fredericton A Halifax Int’l A Charlottetown A St. John's A

BC YT NT NU AB SK MB ON ON QC NB NS PE NL

Climate ID

1018620 2101300 2204100 2402590 3012205 4016560 5023222 6106000 6158733 7016294 8101500 8202250 8300300 8403506

# of Years

46 47 46 46 50 46 46 50 47 42 42 50 40 34

50-year return period value

500-year return period value

MOM

MML

MLM

GLSM

GLSM-V

MOM

MML

MLM

GLSM

GLSM-V

80.5 75.9 67.8 112.9 81.0 97.0 87.6 86.6 95.5 87.9 77.8 103.1 92.9 121.2

84.0 79.8 70.3 118.3 81.3 98.3 85.5 89.0 98.3 91.0 75.4 104.7 94.2 119.3

81.5 76.9 68.9 114.4 81.1 97.8 87.8 87.7 96.6 89.2 76.6 104.3 93.4 119.9

84.2 80.2 70.5 118.7 81.6 98.6 85.6 89.3 98.8 91.5 75.7 105.3 93.8 119.8

84.7 79.8 71.9 116.9 81.0 99.3 86.5 89.7 99.0 92.2 74.2 106.2 94.0 118.3

94.2 89.5 79.0 137.1 94.3 111.7 101.7 100.5 108.6 101.4 90.3 123.0 107.7 141.5

99.9 96.0 83.0 145.8 94.7 113.9 98.2 104.5 113.1 106.4 86.3 125.6 110.0 138.5

95.8 91.3 80.9 139.5 94.5 113.1 102.0 102.4 110.5 103.6 88.3 125.1 108.7 139.3

100.2 96.6 83.4 146.5 95.2 114.4 98.4 105.0 114.0 107.3 86.8 126.6 109.2 139.2

100.8 95.7 85.6 143.3 94.3 115.5 100.0 105.5 114.3 108.2 84.5 128.0 109.5 136.9

using the considered methods; the Eff, Bias and RMSE are calculated for N replications. Before carrying out systematic assessment, numerical analyses were completed for a range of N values and for u ¼ 0 and a ¼1 (i.e., reduced Gumbel variate, Y). It was concluded that the variance of ζ^ (which is required to calculate the Eff shown in Eq. (15)), Bias and RMSE exhibit negligible change for N greater than 50,000. Consequently, the reduced Gumbel variate and N ¼ 50,000 are considered for the systematic assessment. The obtained results are presented in Table 3. From the table, it can be concluded:

Fig. 5. Sites of the selected Canadian meteorological stations.

If Eff is greater than 1, it indicates that the method I is preferred over method II. In other words, the preferred method is the ^ The Bias and RMSE are one associated with a lower variance of ζ. given by, Bias ¼

1 N ^ ∑ ððζ −ζÞ=ζÞ; Ni¼1 i

ð13Þ

and, RMSE ¼

1 N ∑ ððζ^ −ζÞ=ζÞ2 Ni¼1 i

!1=2 ð14Þ

where N is the number of replications and ζ^ i denotes the estimator of ζ for the ith replication. The variance of the estimator of ζ equals (RMSE2−Bias2)ζ2. The performance assessment is carried out using the simple Monte Carlo technique. A set of samples of size n for a replication is generated from the Gumbel model by using, x ¼ u þ að−lnð−lnðwÞÞÞ

ð15Þ

where w is a value of the standard uniformly distributed random variate between 0 and 1. For the set of samples, ζ^ i is calculated

(1) Based on Eff alone (i.e., the variance alone), the MML, which is slightly better than the GLSM, is ranked the best. This is followed by the GLSM-V, MLM and MOM, in descending order of preference. (2) In general and in terms of Bias, the GLSM, GLSM-V and MLM are comparable. All three of these methods are better than the MOM and MML. The MML is most biased among the methods considered. The finding that the MOM is less biased than the MML is in agreement with Lowery and Nash (1970), and that the MLM is less biased than the MOM is agreement with Landwehr et al. (1979). (3) By considering the RMSE alone, the MML and the GLSM are comparable. These two methods are slightly better than the GLSM-V. The MLM is less preferred than the MML, GLSM and GLSM-V, and is more preferred than the MOM. Note that the order of preference according to Eff and RMSE is similar, which is expected as the biases associated with the methods (especially with the GLSM, GLSM-V and MLM) are small. Overall, it is concluded that the GLSM and GLSM-V are the most preferred methods. Also, the analysis leading to the results shown in Table 3 is repeated for u ¼ 65 and a ¼6.5 (i.e., mean of 69 (km/h) and coefficient of variation of 0.12 representing typical annual maximum wind speed statistics). The trends and conclusions that can be drawn from these results are similar to those from Table 3.

4. Application for actual wind records For the application of the discussed distribution fitting methods to historical wind records, we consider the hourly wind data recorded at 14 meteorological stations across Canada. The stations included in the analysis are listed in Table 4 and shown on a map of Canada in Fig. 5. The set of stations is comprised of the capital

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132

cities of each province and territory, as well as the national capital. It is shown in Table 4 that the length of the wind record for these stations is at least 34 years, and for 10 stations it is greater than 45 years. The wind speed records from the Environment Canada (EC) HLY01 database were considered in this study, which consists of hourly aviation (1- or 2-min averages recorded just before the top

127

of the hour) wind speeds or surface weather (10-min average recorded just before the top of the hour) wind speeds dating back to 1953. It is considered that the wind speeds in the HLY01 database are representative of the hourly mean wind speed; this assumption may lead to slightly conservative results (Newark, 1984).

Fig. 6. Fitted Gumbel distributions to the wind speed for the selected meteorological stations.

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132\94\maccounttest15=tx1

Fig. 6. (continued)

Uncertainties exist when considering historical wind data including anemometer siting (i.e. height and surrounding obstacles), anemometer type, and upstream terrain changes over time. Some of these uncertainties and their treatment are discussed by Yip et al. (1995) and Wan et al. (2010). For the current work, anemometer site records were obtained from EC and assessed to determine the quality of wind data at each station. These records included changes in anemometer height, location and instrumentation. The wind speed measurements at each station were standardized to 10 m height using the power law with an exponent of 1/7 (Wan et al., 2010). For some stations, most commonly in the 1953–1970 era, the anemometer was located on the roof of a building. This results in wind speed measurements which cannot be corrected, as they are often greatly affected by the local aerodynamics of the building. In these instances, the affected years of data were removed from the dataset. While the majority of the data used in the current analysis consists of measurements from U2A anemometers, measurements exist from 45B and 78D anemometers. The assessment of the uncertainty due to anemometer type and instrumentation is beyond the scope of this work. Overall exposure corrections based on surrounding terrain conditions were applied to the wind speed data for each location, and the annual maxima were extracted. For the extraction, quality

control techniques were enforced such as a minimum number of wind observations per year (close to 24 observations per day) and verification of extreme maximum wind speeds. No attempt is made to consider uncertainty in climate variability, siting uncertainty, or sampling error of the anemometers. By considering the annual maximum wind speeds as a Gumbel variate, the annual maxima and fitted Gumbel distributions using the various fitting methods are plotted on Gumbel probability paper for each station in Fig. 6. The plots indicate that, for the majority of stations, the Gumbel model seems adequate. In several cases, the MOM and MLM result in similar fits, and the MML, GLSM, and GLSM-V result in similar fits. This grouping trend is not consistent with the results using the simulated samples (as shown in Table 3). This disagreement may be caused by statistical variability or inadequacy of the distribution model. As mentioned previously, on occasion V 2AH instead of VAH is modeled as Gumbel variate. This is based on the consideration that the hourly-mean wind speed can be modeled as Weibull variate, and the extremal distribution of an independent and identically distributed Weibull sequence should be used (Hong, 1994; Cook and Harris, 2004). By considering that V 2AH is a Gumbel variate, probability distribution fitting similar to those shown in Fig. 6 was carried out for the squared wind speed. The fitted distributions are shown in Fig. 7.

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132

Comparison of Figs. 6 and 7 indicate that in some cases the fits of the wind speed are preferred. The estimated 50-year and 500-year return period wind speeds for each location resulting from the different fitting methods are shown in Table 4 by considering VAH as a Gumbel variate. Using the return period wind speeds corresponding to the GLSM as the benchmark, the results indicate that, in general, the MOM and MLM underestimate the return period wind speeds. The underestimation is most severe for the MOM. For all the stations, the

129

maximum underestimation (by the MOM) is 5.5% for the 50-year return period wind speed and 7.5% for the 500-year return period wind speed. The MML underestimates the return period wind speeds only slightly; the relative difference of the estimated return period wind speed by the MML and GLSM is within 2%. The GLSM-V can under- or over-estimate the return period wind speed, depending on the station considered. The relative difference of the estimated return period wind speeds by the MLM and GLSM is approximately 4%. Considering the fact that 34 to 50 years

Fig. 7. Fitted Gumbel distributions to squared wind speed for the considered meteorological stations.

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132\94\maccounttest15=tx1

Fig. 7. (continued)

of data for each station are used for the analysis, the mentioned differences due to distribution fitting alone cannot be neglected. If the V 2AH rather than VAH is modeled as the Gumbel variate, the obtained 50-year and 500-year return period wind speeds (which are calculated based on the return period of V 2AH ) are shown in Table 5. Again, using the return period wind speeds corresponding to the GLSM as the benchmark, the maximum absolute value of the relative difference is 4.5% for the 50-year return period wind speed and 6% for the 500-year return period wind speed. In all cases, the estimated return period values by assigning V 2AH as a Gumbel variate are smaller than those by assigning VAH as a Gumbel variate. This observation is expected and explained in Cook and Harris (2004). The preference and debate on whether the extreme wind speed or its squared value is Gumbel distributed are beyond the scope of this study. Before finalizing this application section, we note that the measured wind speeds were rounded. This rounding error in some cases may affect the estimated return period value of the wind. To quantify this effect, and following Harris (2009), we consider that each measured wind speed has an uncertain rounding error distributed uniformly between [−0.5, 0.5] (km/h). By considering VAH as Gumbel variate and using the GLSM to fit the

historical wind speed including the sampled rounding error, the 50-year and 500-year return period values of VAH were estimated. This process was repeated 10,000 times, the average (Ave.), the minimum (Min.) and the maximum (Max.) values of the return period values of VAH were calculated and are shown in Table 6. The table indicates that the average of the return period values are almost identical to those obtained without considering the rounding error; the difference between the maximum and minimum estimated return period values is less than 1.7 (km/h) for 50-year return period value, and 2.8 (km/h) for 500-year return period value. The standard deviation of the estimates (not shown in the table) is less 0.3 (km/h) for 50-year return period value, and less than 0.4 (km/h) for 500-year return period value. Moreover, to assess the sample size effect on the estimated return period values, we used the GLSM and repeated the analysis that was carried out for Table 3 but taking the fitted Gumbel model obtained by the GLSM to the actual wind speed data as the actual distribution of the extreme wind speed. The obtained results are also shown in Table 6, indicating that in all cases, the estimates are almost unbiased and the (normalized) RMSE is less than 5% for 50-year return period value and 6% for 500-year return period value.

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132

131

Table 5 Estimated return period values of hourly-mean wind speed (km/h) considering that the annual maximum squared wind speed is a Gumbel variate. Location

50-year return period value

Victoria Int’l A Whitehorse A Yellowknife A Iqaluit A Edmonton Int’l A Regina Int’l A Winnipeg Int’l A Ottawa Int’l A Toronto Int’l A Quebec Int’l A Fredericton A Halifax Int’l A Charlottetown A St. John’s A

500-year return period value

MOM

MML

MLM

GLSM

GLSM-V

MOM

MML

MLM

GLSM

GLSM-V

78.6 73.8 66.2 109.4 80.0 95.6 86.6 84.9 94.1 86.2 77.4 100.7 91.2 120.5

80.2 75.4 67.5 109.7 78.7 95.5 83.9 85.7 95.5 87.7 73.9 99.9 91.2 116.0

79.2 74.6 67.1 109.9 79.6 95.9 86.3 85.6 94.9 87.2 75.8 101.1 91.5 118.0

80.3 75.8 67.6 110.3 79.0 95.7 83.8 85.9 95.9 88.1 74.1 100.3 90.9 116.5

80.8 75.8 68.8 109.6 78.5 96.3 84.5 86.4 96.2 88.8 72.9 101.0 91.1 115.3

88.7 83.7 74.5 126.1 90.3 107.0 97.5 95.4 104.4 96.4 87.4 115.1 102.5 136.6

91.0 86.0 76.3 126.6 88.4 106.9 93.6 96.6 106.4 98.6 82.2 114.0 102.5 130.0

89.6 84.8 75.7 126.8 89.7 107.5 97.1 96.5 105.6 97.9 84.9 115.7 102.9 132.9

91.2 86.5 76.5 127.4 88.8 107.2 93.5 97.0 107.0 99.2 82.6 114.5 102.1 130.6

91.8 86.5 78.2 126.4 88.2 108.0 94.6 97.5 107.4 100.2 80.8 115.6 102.4 129.1

Table 6 Effect of rounding error and effect of sample size on the estimated return period values of hourly-mean wind speed (km/h) considering that the annual maximum wind speed is a Gumbel variate (Bias and RMSE are normalized values defined by Eqs. (13) and (14)). Location

Victoria Int’l A Whitehorse A Yellowknife A Iqaluit A Edmonton Int’l A Regina Int’l A Winnipeg Int’l A Ottawa Int’l A Toronto Int’l A Quebec Int’l A Fredericton A Halifax Int’l A Charlottetown A St. John's A

Samples size (i.e., # of Years)

46 47 46 46 50 46 46 50 47 42 42 50 40 34

Effect of rounding error

Effect of sample size

50-year return period value

500-year return period value

50-year return period value

500-year return period value

Ave.

Min.

Max.

Ave.

Min.

Max.

Bias (  100)

RMSE

Bias (  100)

RMSE

84.3 80.3 70.6 118.8 81.6 98.7 85.7 89.4 98.8 91.6 75.8 105.3 94 119.9

83.5 79.6 69.9 118.1 81 98.1 85 88.6 98 90.8 75 104.6 93 119

85 81 71.3 119.4 82.2 99.3 86.2 90.1 99.6 92.3 76.6 106 94.7 120.7

100.4 96.7 83.6 146.6 95.4 114.6 98.5 105.2 114 107.4 87 126.7 109.6 139.4

99.1 95.6 82.4 145.5 94.3 113.6 97.4 103.9 112.6 106.1 85.7 125.6 108.1 137.9

101.6 97.9 84.7 147.5 96.3 115.6 99.4 106.4 115.3 108.6 88.2 127.8 110.7 140.7

−0.01 −0.02 0.01 0.01 0.01 0.02 0.00 −0.02 0.00 0.01 0.01 −0.01 −0.01 −0.01

0.04 0.05 0.04 0.05 0.04 0.04 0.03 0.04 0.03 0.04 0.04 0.04 0.04 0.04

−0.01 −0.03 0.01 0.02 0.01 0.03 0.00 −0.03 0.00 0.01 0.02 −0.01 −0.01 −0.01

0.05 0.06 0.05 0.06 0.05 0.05 0.04 0.05 0.05 0.05 0.05 0.05 0.05 0.06

5. Conclusions The Gumbel distribution is often used to model extreme wind speeds. The most commonly used fitting methods are the MOM, MML and MLM. The use of the generalized least-squares method (GLSM) (i.e., Lieblein BLUE) for the fitting is less common, partly due to the unavailability of the coefficients of the BLUE for sample size n greater than 24. Calculation of the coefficients of the BLUE is carried out for n≤100, although there was no difficulty in evaluating the coefficients for larger values of n. It is shown that the normalized coefficients are well behaved and can be approximated by empirical equations. Results from the simulation analysis indicate that in terms of efficiency, Bias and RMSE: 1) Based on efficiency alone, the MML, which is slightly better than the GLSM, is ranked the best. This is followed by the GLSM-V, MLM and MOM, in descending order of preference. 2) In general and in terms of Bias, the GLSM, GLSM-V and MLM are comparable. All three of these methods are better than the MOM and MML. This preference is followed by MOM and MML in descending order. 3) By considering the RMSE alone, the MML and the GLSM are comparable. They are slightly better than the GLSM-V. This preference is followed by the MLM and then by the MOM.

Overall, it can be considered that the GLSM and GLSM-V are the most preferred methods. As the evaluation of the coefficients of the BLUE in the GLSM is more involved than that in the GLSM-V, the latter could be used as a substitute for the former. Application of the fitting methods for annual maximum wind speeds at several Canadian sites is presented. It was observed that the fitted model by the MOM and MLM are similar, while those by the MML, GLSM and GLSM-V are similar if the annual maximum wind speed is considered to be a Gumbel variate. In comparison with the results obtained by using the GLSM, the MOM and MLM underestimate the return period values; the MML also underestimates the return period values, but very slightly; and the GLSM-V can under- or over-estimate the return period value depending on the station considered. These observations for the selected meteorological stations may not be generalized for other stations because of the statistical variability.

Acknowledgments The authors are grateful to R. J. Morris from Environment Canada for providing the historical wind data for the meteorological stations included in this work, as well as providing helpful advice on its interpretation. Useful discussions with D. A. Gatey

H.P. Hong et al. / J. Wind Eng. Ind. Aerodyn. 119 (2013) 121–132\94\maccounttest15=tx1

regarding the processing of historical wind records are greatly appreciated. Financial support received from Natural Science and Engineering Research Council of Canada and the University of Western Ontario is much appreciated. We thank the anonymous reviewers for their constructive comments and suggestions. Appendix. Moments of order statistics for the Gumbel model The moments of order statistics that need to be evaluated are, Z ∞ n! −x −x Eðyki:n Þ ¼ xk e−x−ie ð1−e−e Þn−i dx ðA1Þ ði−1Þ!ðn−iÞ! −∞ and Eðyi:n yj:n Þ ¼

n! ði−1Þ!ðj−i−1Þ!ðn−jÞ! Z ∞Z ∞ −y −x −y −x −y  xye−y−e −x−ie ðe−e −e−e Þj−i−1 ð1−e−e Þn−j dxdy −∞

−∞

ðA2Þ Explicit equations to evaluate derived by Lieblein (1953), are Eðyki:n Þ ¼

Eðyki:n Þ

and Eðyi:n yj:n Þ, which are

n−i n! ðn−iÞ! ∑ ð−1Þr g ði þ rÞ ði−1Þ!ðn−iÞ! r ¼ 0 r!ðn−i−rÞ! k

ðA3Þ

and n! ði−1Þ!ðj−i−1Þ!ðn−jÞ! j−i−1 n−j ðj−i−1Þ! ðn−jÞ! hði þ r; j−i−r þ sÞ  ∑ ∑ ð−1Þrþs r!ðn−i−1−rÞ! s!ðn−j−sÞ! r¼0s¼0

Eðyi:n yj:n Þ ¼

ðA4Þ

where g k ðcÞ ¼ ð−1Þk dtdk ðΓðtÞc−t Þjt ¼ 1 , 2tuhðt; uÞ ¼ ðu−tÞg 2 ðt þ uÞ þ t 2 ðg 1 ðtÞÞ2 þ 2Lð1 þ t=uÞ−ðlnðt=uÞÞ2 −π 2 =6;

ðA5Þ

and Γ() is the gamma function, and Lð1 þ t=uÞ is the Spence's integral (or the dilogarithm integral) (Abramowitz and Stegun, 1972). Our experience indicates that the use of the explicit equations for numerical evaluation is suitable only for n less than about 20 because of the alternative signs in Eqs. (A3) and (A4). Balakrishnan and Chan (1992a) used numerical integration for their evaluation and provided the covariance matrix for n up to 30. For this study, we use the Gauss–Kronrod quadrature implemented in MATLAB (Shampine, 2008) to solve Eqs. (A(1) and A2). We verified the numerical values with those calculated by using Eqs. (A3)–(A5) for n up to 15 and, with those reported by Balakrishnan and Chan (1992a) for n up to 30. The verification is also carried out using the following identity (David, 1981), n

∑ Eðyki:n Þ ¼ nEðY k Þ; for k ¼ 1 and 2

i¼1

ðA6Þ

and n

n

∑ ∑ βi;j:n ¼ nV arðYÞ

i¼1j¼1

where βi;j:n ¼ EðY i:n Y j:n Þ−αi:n αj:n . Values of cu and ca are estimated by using Eq. (10).

ðA7Þ

References Abramowitz, M., Stegun, I.A. (Eds.), 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing. Dover, New York. Balakrishnan, N., Chan, P.S., 1992a. Order statistics from extreme value distribution, I tables of means, variances and covariances. Communication in Statistics— Simulation and Computation 21 (4), 1199–1217. Balakrishnan, N., Chan, P.S., 1992b. Order statistics from extreme value distribution, II best linear unbiased estimates and some other uses. Communication in Statistics—Simulation and Computation 21 (4), 1219–1246. Castillo, E., 1988. Extreme Value Theory in Engineering. Academic Press Inc., New York. Cook, N.J., 1985. The Designer's Guide to Wind Loading of Building Structures Part 1. Butterworths, London. Cook, N., Harris, R.I., 2004. Exact and general FT1 penultimate distributions of extreme wind speeds drawn from tail-equivalent Weibull parents. Structural Safety 26, 391–420. David, H.A., 1981. Order Statistics, second ed. Wiley, New York. Draper, N.R., Smith, H., 1998. Applied Regression Analysis, third ed. Wiley, New York. Frank, H., 2001. Extreme Winds over Denmark from the NCEP/NCAR Reanalysis. Risø National Laboratory, Technical Report Risø-R-1238(EN). Harris, R.I., 1996. Gumbel re-visited—a new look at extreme value statistics applied to wind speeds. Journal of Wind Engineering and Industrial Aerodynamics 59, 1–22. Harris, R.I., 2009. XIMIS, a penultimate extreme value method suitable for all types of wind climate. Journal of Wind Engineering and Industrial Aerodynamics 97, 271–286. Hong, H.P., 1994. A note on extremal analysis. Structural Safety 13 (4), 227–233. Hosking, J., 1990. L-moments: snalysis and estimation of distributions using linear combinations of order statistics. Journal of the Royal Statistical Society. Series B. Statistical Methodology 52 (1), 105–124. Hosking, J., Wallis, J., Wood, E., 1985. Estimation of the generalized extreme-value distribution by the method of probability-weighted moments. Technometrics 27 (3), 251–261. Landwehr, J.M., Matalas, N.C., Wallis, J.R., 1979. Probability weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resources Research 15, 1055–1064. Leadbetter, M.R., Lindgren, G., Rootzen, H., 1983. Extremes and Related Properties of Random Sequences and Processes. Springer-Verlag Inc., New York. Lieblein, J., 1953. On the exact evaluation of the variances and covariances of order statistics in samples from the extreme-value distribution. The Annals of Mathematical Statistics 24 (2), 282–287. Lieblein, J., 1974. Efficient Methods of Extreme-value Methodology. National Bureau of Standards, Washington, Report NBSIR 74-602. Lloyd, E.H., 1952. Least-squares estimation of location and scale parameters using order statistics. Biometrika 39, 88–95. Lowery, M.D., Nash, J.E., 1970. A comparison of methods of fitting the double exponential distribution. Journal of Hydrology 10, 259–275. Newark, M.J. 1984. Canadian wind data. In: Proceedings of the Fourth Canadian Workshop on Wind Engineering, Canadian Wind Engineering Association, Toronto, Canada. Peterka, J.A., Shahid, S., 1998. Design gust wind speeds in the United States. Journal of Structural Engineering 124, 207–214. Phien, H.N., 1987. A review of methods of parameter estimation for the extreme value type-1 distribution. Journal of Hydrology 90, 251–268. Sacré, C., 2002. Extreme wind speed in France: the ’99 storms and their consequences. Journal of Wind Engineering and Industrial Aerodynamics 90, 1163–1171. Shampine, L.F., 2008. Vectorized Adaptive Quadrature in MATLAB. Journal of Computational and Applied Mathematics 211, 131–140. Wan, H., Wang, X.L., Swail, V.R., 2010. Homogenization and trend analysis of Canadian near-surface wind speeds. Journal of Climate 23, 1209–1225. Yip, T., Auld, H., Dnes, W., 1995. Recommendations for updating the 1995 National building code of Canada wind pressures. In: Proceedings of the Ninth International Conference on Wind Engineering. IAWE, New Delhi, India.