Statistics of extreme values of air quality—a simulation study

Aawsphrlc

Enviromunr Vol. 19, No. IO. pp. 1713-1721.

00044981

1985

85 13 00 l 0.00

Pergamon Res

Printed in Great Bnuia.

Ltd.

STATISTICS OF EXTREME VALUES OF AIR QUALITY-A SIMULATION STUDY DAVID P. CHOCK Environmental Science Department, General Motors Research Laboratories, Warren, MI 48090, U.S.A. (First received 3 October 1984 and receioedjor publication 18 March 1985) Abstract-Some of the existing National Ambient Air Quality Standards require the use of an extreme observed concentration in a year to determine compliance. Since observed extreme values tend to be not reliable, different statistical approaches for determining the extreme values have been used or suggested. However, none of these approaches properly take into account the effects of an underlying trend and the serial correlation of the air quality time series. By means of a time series simulation, these effects can be considered concurrently in estimating the extreme values. This paper reports the results of such a simulation for determining the statistics of the seven highest values (rank m = l-7, m = 1 representing the highest value) using actual air quality data that contain both trends and autocorrelations. The result of this simulation shows that for a high-pollution season of 122 days, the commonly used asymptotic distributions overestimate the maximum (m = 1)values and underestimate their uncertainties. As one moves from m = 1 to m = 7, the over- and underestimations by the asymptotic distributions worsen (compared to the simulation result). These findings in logarithmic space are further enhanced when they are converted back to concentration space. The simulation using the oxidant data for Azusa, California further shows that the uncertainties associated with the estimates of the extreme values are typically 20 y0 of the values form = 1and 10 y0of the values form = 7. Compared to the observed data, which is a single series for each year, the result based on the popular lognormal distribution consistently overpredicts the extreme values, by about 40% for the maximum values and about 20% for the seventh highest values. Our results illustrate the difficulty of estimating the extreme values of air quality time series with accuracy and confidence. However, the accuracy and confidence of the estimates improve as the rank moves away from the extreme. This result calls for the need for using a less extreme value in setting a sensible air quality standard. Of course, such a standard can be set without sacrificing its stringency.

1. INTRODUCIION

upper extremes of air quality data have generally been used both to ascertain the attainment of air quality standards and to implement a specific emission control strategy once a standard is violated. In the U.S., some of the National Ambient Air Quality Standards have been set or proposed such that a standard concentration for a specific pollutant can be exceeded no more than once per year (on the average). If the expected number of exceedances is greater than one The

per year, then the emissions which lead to the observed air quality should be adjusted so that the second-

highest concentration of the year does not exceed the standard concentration in the future. It is generally realized that the extreme values of air quality data may well be statistical outliers and may contain much uncertainty. Therefore, it is precarious to accept the observed concentrations at or near the upper extreme at face value and implement control strategies based on them. To determine the upperextreme values, what has often been done is to take one year’s worth of data and fit a statistical distribution to them. The estimated extreme concentrations can then be determined according to their corresponding cumulative probabilities in the resulting distribution. This approach may mitigate the outlier problem. But, since

the data do not necessarily fit one statistical distribution unambiguously, and since different distributions have very different behavior near the extreme tails, estimates of the extreme values can have substantial uncertainty due to the possibility of fitting different distributions to the data. To reduce this uncertainty, a statistical distribution is often fitted to only the top 20-50 % of the data so that the extreme values in effect become ‘less extreme’ in the fitted distribution(s). The above approaches contain many unwarranted assumptions. First of all, fitting a statistical distribution to a data set presumes that each of the data points is an independent, identically distributed (i.i.d.) random variate. But in actual fact, it is well known that air quality data are generally correlated (see, e.g. Chock et ol., 1975; Hirtzel and Quon, 1981; Chock, 1982). Moreover, they also exhibit seasonal trends. These two properties negate the i.i.d. assumption. One can alleviate the seasonal-trend problem by considering the data from the high-pollution season only. Even so, deterministic trends are often evident. As an example, Fig. 1 shows the time series of the logarithms of the daily maximum l-h oxidant concentrations (in pphm) in Azusa, California from 1 June to 30 September 1970. If we simply ignore the deterministic trend, which is approximated by a cubic polynomial (to be discussed later) and displayed as a solid line in Fig. 1, the

1713

DAVID P. CHOCK

.

.

* .

.

. .

.

.

. .

.

.

*.

.

.

.

.

.”

-.

.

.

.

Fig. 1. Time series of the logarithms of the daily maximum l-h oxidant concentrations (in pphm) in Azusa. California from I June lo 30 September 1970. The solid line is the polynomial trend whose coefkients arc displayed in fable I.

cumulative probability plot (on a normal probability scale) of the data after removing the series mean is shown in Fig. 2(a) by closed circles. However, if the trend is removed, the residuals will acquire a different distribution pattern, shown in open circles in Fig. 2(b). Since the detrended data points are correlated (the autocorrelation for lag of one day, p,, is 0.446 in this case), one can ask what the distribution would be like if the autocorrelations were removed. A simulated example, which assumes a normal distribution with an

appropriate variance (to be discussed later), is shown as crosses in Fig. 2(c). While it is obvious that for a specific data set, the distributions for trended and detrended series are likely to be different, it is not obvious how the cumulative probability distribution of data from a detrended, correlated time series should differ from that of a detrended but uncorrelated series. Strictly speaking, data of a correlated series should be characterized by a multivariate distribution (see Chock, 1984) and the cumulative probability plot in

1.2 c 09

-

fc

0.6

-

z

0.3

-

z k .c-

#

n/

o-

8t” .0.3

-

;

-0.6

-

-0.9

-

-1.2

-

I 0.01

I

I

I

lo.11 0.5 I 2 0.050.2 1

1 5

10

I 20

I

I

30

I50 I70 40 60

Probability

I

I 60

I

I 90

I 95

98

I

I 99.d 99

J

99.99

99.9

in %

Fig. Z(a). Cumularivc probability plot of the time series (trended and correlakd) described in Fig. 1, but with the seriesmean removed. The solid line, shown in Figs 2(b) and 2(c) as well, is drawn to facilitate comparisons among Figs 2(a)-(c).

Statistics of extreme values of air quality

::

0

.E (j

-0.3

s ’ F

-0.6

0 0.050.2

1

5

20

40

60

60

95

99

99.9

Probability in %

Fig. 2(b). Cumulative probability plot of the time series described in Fig. 1, but with the trend removed.

0.050.2

1

5

60

80

Probability

in %

20

40

95

99

99.9

Fig. 2(c). Cumulative probability plot ofa detrended, uncorrelated time series with zero mean and the same variance as the series described in Fig. 2(b).

Fig. 2(b) is not a marginal variate probability

distribution of the multidistribution. Hence, the statistical

meaning of Fig. 2(b) is not clear, and comparing Figs 2(b) and (c) can lead to misleading conclusions. Second, by collecting only the upper portion of the data the roles of the trend and autocorrelation on the statistical distribution can no longer be ascertained. Even a proper time series analysis of the data becomes impossible because the data set is no longer qually spaced in time. Third, determining the extreme values by reading off

the values of the corresponding quantiles in the fitted distribution ignores the statistical property of extreme values. It is not even obvious what the uncertainties of the estimated extreme values are unless the distributions of the extreme values are considered. Applications of extreme-value statistics to air quality data have been carried out by many investigators (e.g. Roberts, 1979; Johnson and Symons, 1980). But the effect of trend and autocorrelation on extreme value distributions is not well known. There have been studies dealing separately with the extreme values of a

1716

DAVID P. CHOCK

series (e.g. Horowitz and Barakat. 1979; Horowitz. 1980; Chock, 1984) and with the extreme values of a correlated series without a trend (e.g. Hirtzel et 01.. 1982; Chock, 1984, and references in both). Horowitz and Barakat (1979) derived the asymptotic distribution of the maximum value, assuming that the effect of autocorrelation is negligible in the asymptotic limit (i.e. the number of observations is large, as in a sufficiently long series). Chock (1984) derived the exact distributions of the first four highest values of a first-order Markov normal process. The result shows that an increase in autocorrelation will decrease the extreme-value estimates and increase the uncertainties of the estimates, but that the effect of autocorrelation becomes small as the number of observations increases. Chock (1984) also shows the asymptotic distributions of the first m highest values of an uncorrelated series whose mean and variance are functions of time. There has been no study on the effect of both trend and autocorrelation on extreme-value distributions. The present work represents an attempt in that direction. Chief among the questions to be answered are whether the asymptotic (large number of observations) extreme-value distributions are appropriate to provide extreme-value estimates, and whether there are better alternatives to the present criteria used in the air quality standards. We shall address the first question in detail Our approach is to use time-series simulation in order to reduce the constraints placed on the form of the time series. Parameters from actual air quality data were used to generate many series. The generated series were then used to investigate the property of their extreme values. Section 2 describes the method of analysis. The results and discussion are presented in Sections 3 and 4, respectively. In Section 4, the question of a more robust air quality standard is briefly discussed. trended

2. *METHOD OF

ASALYSIS

To determine the extreme values of a correlated series with a deterministic trend, the series should first be detrended. An appropriate model, i.e. Box and Jenkins (1976). can then be determined for the detrended series. The model is then used to generate a large number of series of the same length. The deterministic trend is added back to each of these series. Values of equal ranks in the upper extreme of the resulting series are then averaged and their variances determined. We applied this approach to the logarithms of the daily maximum l-h oxidant concentrations (in pphm) in Azusa, California from 1 June to 30September 1969, until 1978. Within this four-month period of each year, there are 122 consecutive observations. Only the high-pollution season (June through September for ozone and oxidant) was considered for three reasons: (1) to minimize mixing data of different characteristics, (2) to minimize the influence of lower

concentrations on the distribution used to estimate the concentrations at the upper extreme and (3) to mitigate the ‘extremeness’ of the extreme values so that their estimates are somewhat less dependent on the behavior of the parent distribution. Logarithmic data were used for two reasons. One was that a lognormal distribution is commonly used or assumed to fit the air quality data, even though fitting it toa truly detrended, uncorrelated air quality data set has never been done. The other was that the white noise generated in the simulation process is normally distributed so that taking the antilogarithms of the white noise will lead to a lognormal series. The data do not show any clear weekly trend for any year. However, deterministic trends which are due primarily to seasonal variation of the pollutant are evident. An example of the trend was already shown in Fig. 1. To remove the trend, we chose a polynomial fitting approach rather than the time-differencing approach for the following reasons. First, the trend was not stochastic (Box and Jenkins, 1976) within the time period considered. Second, we were interested in generating series which had identical trend behavior. Forecasting was not our intention here. Our fitted , trend has the form x u, rJ. We further scaled the time “iI

factor so that a period of 365 days corresponds to 2 n radians. We also set the time origin at the sixty-first observation (31 July, essentially the midpoint of our time period). So, t = (day of observation in the series61) n/l 83, day one being 1 June. We assumed a cubic polynomial (r = 3). This was found to be adequate. A higher-order polynomial tended to introduce spurious oscillations; a lower one did not fit the trends well in a few instances. Our hope in detrending was to remove only the seasonal cycle while leaving smaller-scale fluctuations alone as much as possible. Since we were using the complete time series to determine the trend rather than using a small portion of it at a time as in the moving-average approach, only the long-period oscillations were filtered, and furthermore, the introduction of extraneous oscillations (Slutzky-Yule effect (see Kendall, 1976; Mihram, 1972)) was also minimized. The estimated coefficients of the trends are shown in Table 1. After the trend from each series was removed, a Box-Jenkins (1976) model was identified for the resulting series. This model has the general form l#$@)O,(B’)VdVPZ,

= 61&B)8Q(B3n,+efJ,

(1)

where Z, and o, are the variates of the time series and of the white noise, respectively, at time t. go is the deterministic trend. B and B’ are the backward shift operators such that B’X, = X,_, and (B’)‘= X,_,,, s being the seasonal period. 4,(B) = P mP(B’)

=

1-

1 -

f

4,

B'

and

i=l

c @,B” are the autoregressive oper,= I ator of ord:r pand seasonal autoregressive operator of

Statistics of extreme values of air quality Table 1. Estimated coefficients of the trend a0

1969 1970 1971 1972 1973 1974 1975

3.328 (0.0617)t 3.430 (0.0549) 3.274 (0.0592) 2.845 (0.0681) 2.996 (0.0630) 2.966 (0.0507) 2.842 (0.0649)

1976 1977 1978

2.783 (0.0706) 2.972 (0.0614) 3.011 (0.0597)

0.437 (0.1701) 0.346 (0.1514) -0.002* (0.1632) -0.451 (0.1876) -0.605 (0.1735) 0.238 (0.1396) 0.153 (0.1789) -0.299 W7) (0.1692) - 0.482 (0.1644)

- 0.696 (0.1260) -0.591 (0.1121) -0.821 (0.1209) -0.173 (0.1389) -0.393 (0.1285) - 0.072 (0.1034) -0.195 (0.1325) -0.479 (0.1442) -0.522 (0.1253) -0.261 (0.1218)

0.131. (0.2371) - 0.329 (0.2109) 0.122. (0.2274) 0.407 (0.2614) 0.870 (0.2418) -0.151 (0.1946) 0.039* (0.2493) 0.027’ (0.2714) -0.485 (0.2358) 0.348 (0.2291)

*Value which is insignificant at the 5 “/, level. tstandard error of the estimate.

order P, respectively. 0,(B) = 1 -

2 BiBiand O,(Bs) i=*

5i

1 - i$t OiB” are the moving-average operator of

order q and the seasonal moving-average operator of order Q, respectively. Vd = (1 - B)d and VP = (1 - Bs)D are the difference operator of order d and the seasonal difference operator of order D, respectively. The general notation for the model is (p,d,q) x (P,D,Q), or just (p,d,q) where applicable. For our detrended series, B. = 0. For d and D = 0, Z, can be replaced by Z, = Z, - p where p is a constant, equal to zero in our case. To help identify a model for a series, autocorrelations, inverse auto-correlations (Cleveland, 1972) and partial autocorrelations were plotted. Generally for each series, two or more models were tested and their parameters were estimated using nonlinear leastsquares iterations. A model was selected based on the chi-square test on the residual autocorrelations. If the residuals a, are uncorrelated white noise, then their

autocorrelations should be normally distributed. (See Mihram. 1972, p. 474.) A brief description of how the model was selected for each year is presented in Appendix A. The selected models, together with their parameters are shown in Table 2. This table also contains ef and af for each year. They are the variances of the detrended series and of the residuals. The latter is modeldependent and is assumed to be the variance of the white noise in the time-series simulation. The autocorrelations with the lag of 1 day, pi, of the detrended series are also included in Table 2. Once a model was established for the data of a specific year, it was used to generate many series which simulated the original data. The white noise in each series was generated (IMSL, 1980) with an appropriate mean and variance. The mean was taken as zero, the white noise variance was supplied by af shown in Table 2. The white noise variance is actually related to the series variance 0: (Box and Jenkins, 1976). For an autoregressive series of order p, u: = (1 - pt +t - . . - ppq$,)u: where p, is the autocorrelation at lag i. For a moving-average series of order q, u: = (1 + 0: + . . . + 0:) ai. The relation is more complicated for an integrated autoregressive-moving-average series. The seasonal model can, of course, be written in terms of a high-order (p,d,q) model, with many intervening parameters being set explicitly equal to zero. The whitenoise generator was first tested to see if its generated series could be safely assumed to be white-noise. The

significance level from the chi-square test on the autocorrelation of the first 24 lags of each of many generated series tested was close to 1. Furthermore, the input white-noise variance was recovered to within 0.5 % on the average. It was concluded that the whitenoise generating capability of our simulator was satisfactory. In the time series simulation, initial values are required for series with p > 0. We used zeros as the starting values when p is small, say, < 2. For a model with a large p. or more generally, for a seasonal model where (Ps + p) starting values were needed, zeros were again supplied as starting values, but the series was allowed to evolve for at least twice the period of (Ps + p) for the transient effect of the starting values to

Table 2. Time series models and their parameters Parameters 1969 1970 1971 1972 1973 1974

1975 1976 1977 1978

1717

0, = -0.699, u: = 0.200, u: = 0.125, p, = 0.581 4, = 0.413, UJ, = -0.275, u2 = 0.158, 6’ = 0.121, p, = 0.446 +r = 0.555, @, = -0.251, CT;= 0.184, uf = 0.126, p, = 0.537 4, = 0.404, c: = 0.243. CT:= 0.207, P, = 0.401 4, = 0.470, $J~= -0.111, 4, = 0.178, 44 = -0.282. u; = 0.208, u: = 0.162, p, = 0.426 4, =0.441,u:=0.134,u:=0.110.p, =0.435 4r = 0.769, I& = -0.231. a: = 0.221, u: = 0.131, p, = 0.624 t’, = -0.753, 0* = -0.404, u: = 0.262, u: = 0.154, p, = 0.618 0, = -0.579, u; = 0.198, uj = 0.155, p, = 0.387 4, = 0.590, 0, = -0.261, CT’= 0.186, at = 0.110, 0, = 0.626

1718

DAVID P. CHOCK

decay before the simulated series was registered. For each year, a total of 100 series were generated, each with 122 observations. For checking fiurposes, in a few instances, the model-identification routine was applied to two or three members of the 100 series. The resulting parameters were qualitatively similar to the original parameters used for the simulation. For series containing seasonal parameters, the oscillatory behavior of the autocorrelation function of the simulated series tended to be more prominent than that in the original series. Notice, however, that the seasonal parameters in Table 2 are all small. Their impacts on the extreme values are likely to be limited. The average variance of the 100 series is within -+ 8 Y<(more typically within & 4 y(L)of the variance of the original detrended series. For the autoregressive models which required initial input values, the average variances are within f 4?/, of their original series variances. ‘Thus, our handling of the initial input values should be considered adequate. Another factor justifying our procedure for initial input values is that pi’s are not too large (see Table 2)and the series are sufficiently long. A persistence problem may occur when pi is large ( 2 0.8, say) and/or when the simulated series is too short, resulting in an underestimated variance for the simulated series. For each year, the highest values from all the series were grouped together and their mean and variance were determined. We chose the mean rather than the mode since the latter was somewhat more cumbersome to determine. Their values approach each other as one moves away from the most extreme values. The second-highest, third-highest, , seventh-highest values from all the series were similarly grouped and their means and variances were determined. The deterministic trend for each year was then added to each of the 100 series of the corresponding year, and the above procedure to determine the means and variances of the extreme values was repeated. In so doing, we can study the extreme values of both trended and detrended series. The results will be discussed in the next section. In addition to the simulation described above, we determined by simulation the same extreme values and their variances but this time ignoring the effect of autocorrelation (i.e. setting all the time-series parameters to zero). For comparison purposes, we also estimated the extreme values using the asymptotic extreme-value distributions for trended but uncorrelated series (Chock. 1984). To do this, numerical integrations of the following type were performed.

where/(x) = .Tfor the mean andf(x) = (II - .&,,)’ for the variance. The subscripts m and n designate the mthhighest value and n observations (n = 122 in our case). For the trended but uncorrelated series, the cumulative asymptotic distribution for the mth highest value in n

observations is (see Appendix of Chock. 1984). (3) where mI)

1

(4)

with 1, (Zlogn)’ *; u,, z - (1/3)((1;2)(log4n + log log n) + log m); 0. standard deviation of the detrended series; P,, the trend value at time i. 3. RESULTS

Before describing the results, it is worthwhile to estimate the uncertainties associated with the means and variances of the highest and the seventh-highest values obtained from 100 (= N) series. From Table 3. which will be discussed below, the estimated variances s* for the samples containing the 100 highest values and the samples containing the seventh-highest values are about 0.04 and 0.01, respectively, based on the data from the detrended series. The corresponding means are typically 1.10 and 0.70, respectively. Therefore, the uncertainties for the sample means and variances are - si JN and 5 s2(2/(N - 1))“’ (strictly speaking, the latter is true only for samples from a normal distribution). or 0.02 (- 2 “/ of sample means) and 0.006 ( - 14Y4 of sample variance), respectively, for the sample containing the highest values, and 0.01 ( - 1.5 Is/,of sample mean) and 0.001 ( - 14 % of sample variance), respectively, for the sample containing the seventh-highest values. These estimates provide us with the typical magnitude of uncertainties of our estimated values. Tables 3 and 4 show the averages and standard deviations of the seven highest values obtained from our simulation of the trended and detrended correlated series for each year. Several points are worth noting. (1) If we determine the maximum of the trend from Table 1,and add it to the average maximum value of the detrended series, we find that the result exceeds the average maximum of the trended series. For example, for 1970, the trend-line maximum is 3.373, the average maximum of the detrended series is 1.024. The sum, 4.397, exceeds the mean maximum of the trended series, which is 4.361. This is expected because the occurrence of the maximum of the detrended series and of the maximum of the trend do not necessarily coincide. (2) Because the trend is not constant, the variance of the trended series is greater than that of the detrended series. Accordingly, the uncertainty of the mean mth-highest value of the trended series tends to be greater than that of the detrended series. (3) Compared to actual observations (for m up to 4, see Table 3 of Chock, 1984). the estimated means of the extreme values are consistently higher than the observed values of the same ranks. The excesses in the

1719

Statistics of extreme values of air quality Table 3. Averages and standard deviations (in parentheses) of the seven highest values from simulation of detrended* correlated series Year

Max

2nd

3rd

4th

5th

6th

7th

1969 1970

I .082 (0.169) 1.024

0.953 (0.139) 0.902

0.874 (0.114) 0.820

0.820 (0.102) 0.762

0.779 (0.099) 0.722

0.741 ‘;;;)

0.706 (2;)

1971

(0.158) 1.105

(0.126) 0.972

(0.113) 0.896

(0.095) 0.830

(0.087) 0.788

(0.081) 0.751

(%

(0.177) 1.251 (0.197) 1.174 (0.184)

(0.139) 1.102 (0.167) 1.051 (0.155)

(0.124) 1.009 (0.146) 0.952 (0.121)

(0.116) 0.940 (0.134) 0.887 (0.110)

(0.111) 0.887 (0.120) 0.838 (0.097)

(0.107) 0.844 (0.112) 0.795

(0.101) 0.807 (0.105) 0.756 (0.089)

(0.149) 0.931 1.176 (0.209) 1.277 (0.213) 1.141 (0.170) 1.055 (0.185)

(0.125) 0.823 1.048 (0.160) 1.140 (0.184) 1.003 (0.148) 0.929 (0.132)

(0.109) 0.754 0.958 (0.133) 1.041 (0.148) 0.923 (0.120) 0.853 (0.116)

(0.101) 0.703 0.899 (0.131) 0.974 (0.140) 0.858 (0.108) 0.790 (0.110)

(0.093) 0.664 0.846 (0.130) 0.916 (0.138) 0.816 (0.103) 0.749 (0.108)

(0:087) (;zj’ 0.806 (0.120) 0.873 (0.135) 0.776 (0.096) 0.714 (0.102)

(k%3;) 0.773 (0.116) 0.835 (0.129) 0.744 (0.092) 0.679 (0.098)

1972 1973

1974 1975 1976 1977 1978

*The detrended series have a mean of zero.

Table 4. Averages and standard deviations (in parentheses) of the seven highest values from simulation of trended* correlated series Year

Max

2nd

3rd

4th

5th

6th

7th

1969

4.346 (0.168) 4.361 (0.174) 4.202 (0.193) 4.085 (0.200) 4.112 (0.220) 3.903 (0.161) 3.985 (0.218) 3.958 (0.238) 3.996 (0.199) 4.048 (0.197)

4.209 (0.131) 4.220 (0.148) 4.076 (0.161) 3.925 (0.161) 3.963 (0.162) 3.787 (0.139) 3.852 (0.175) 3.824 (0.184) 3.856 (0.161) 3.907 (0.138)

4.119 (0.124) 4.138 (0.127) 3.995 (0.154) 3.825 (0.138) 3.870 (0.138) 3.713 (0.114) 3.761 (0.137) 3.727 (0.162) ‘3.769’ (0.134) 3.825 (0.116)

4.055 (0.110) 4.084 (0.107) 3.936 (0.141) 3.761 (0.136) 3.807 (0.114) 3.667 (y$

4.007 (0.099) 4.026 (0.107) 3.885 (0.128) 3.698 (0.123) 3.752 (0.097) 3.629

3.971 (0.100) 3.977 (0.096) 3.836 (0.117) 3.653 (0.114) 3.713 (0.092) 3.598 (0.095) 3.608 (0.126) 3.544 (0.144) 3.615 (0.107) 3.684 (0.104)

3.932 (0.093) 3.940 (0.089) 3.795 (0.110) 3.610 (0.107) 3.672 (0.089) 3.567 (0.088) 3.574 (0.125) 3.503 (0.137) 3.578 (0.101) 3.649 (0.102)

1970 1971 1972 1973 1974 1975 1976 1977 1978

(0:134) 3.651 (0.160) 3.711 (0.126) 3.774 (0.116)

(EE’ (0.132) 3.592 (0.152) 3.661 (0.126) 3.723 (0.110)

ODetrended series plus trend. logarithmic scale are typically 0.2-0.5 for the maxima, 0.1-0.3 for the seventh-highest values, and somewhere in between for the rest. While this is not a proper comparison since the observed values come from only one series per year, the persistence of the pattern is cause for concern. It is known that the two-parameter lognormal probability density function has a long tail (compared to, say, the Weibull distribution of the same mean and variance) and thus tends to overpredict the mean extreme concentrations. Table 5 shows the results of simulation of the trended, uncorrelated series. In these series, except for

the trend, each observation is treated as an independent observation with a variance identical to that of the corresponding detrended, correlated series. Comparing the simulated results of Tables 4 and 5, we see that the series with autocorrelation have smaller maximum values than the series without autocorrelation, but the differencesare no more than 2 y0of the mean values. The differences gradually diminish and the pattern becomes random as we move away from the maximum values. The standard deviations of the correlated series are generally greater than those of the uncorrelated series, as anticipated. The relative dif-

1720

DAVID P. CHOCK

Table 5. Averages and standard

deviations (m parenthcscs) of the seven highest values from simulation of trended, uncorrclatcd s&es

YCU

Max

2nd

3rd

4th

5th

6th

7th

1969

4.430 (0.225) 4.373 (0.202) 4.253 (0.200) 4.124 (0.219) 4.147 (0.193) 3.928 (0.169) 4.035 (0.210) 4.050 (0.225) 4.032 (0.214) 4.132 (0.194)

4.279 (0.161) 4.220 (0.144) 4.082 (0.151) 3.942 (0.151) 3.966 (0.142) 3.804 (0.131) 3.874 (0.167) 3.850 (0.165) 3.852 (0.154) 3.969 (0.133)

4.176 (0.128) 4.136 (0.116) 3.988 (0.123) 3.828 (0.120) 3.861 (0.116) 3.721 (0.102) 3.772 (0.128) 3.732 (0.131) 3.770 (0.123) 3.862 (0.114)

4.108 (0.123) 4.06: (0.103) 3.911 (0.104) 3.752 (0.109) 3.790 (0.101) 3.668 (0.093) 3.709 (0.119) 3.650 (0.120) 3.689 (0.112) 3.793 (0.093)

4.05 I (0.108) 4.02 1 (0.096) 3.861 (0.101) 3.692 (0.094) 3.733 (0.090) 3.619 (0.084) 3.644 (0.111) 3.587 (0.101) 3.633 (0.103) 3.744 (0.084)

4.013 (0.108) 3.986 (0.096) 3.823 (0.100) 3.644 (0.085) 3.691 (0.081) 3.588 (0.079) 3.598 (0.098) 3.535 (0.095) 3.594 (0.098) 3.702 (0.077)

3.97 1 (0.101) 3.941 (0.092) 3.783 (0.088) 3.606 (0.078) 3.651 (0.073) 3.556 (0.077) 3.563 (0.099) 3.494 (0.088) 3.559 (0.093) 3.667 (0.069)

1970 1971 1912 1973 1974 1975 1976 1977 1978

ferences increase as the rank becomes less extreme, reaching 25 % or more for the seventh-highest values. But this pattern is not persistent, especially for the standard deviations of the average maximum values. Several factors contribute to the lack of consistency. Foremost among them are the greater uncertainties associated with our standarddeviation estimates (compared to those of the means) by simulation and the fact that uncertainties of the estimates increase as the rank of interest becomes more extreme. In addition, inclusion of the deterministic trend also distorts the uncertainties of the estimates. The results based on the asymptotic approximation of extreme values of trended, uncorrelated series are shown in Table 6. Here again, for each year, the Table 6. Averages and standard values based on asymptotic

prior to the introduction of the series variance of the corresponding detrended, correlated series. Comparing Table 6 with Tables 4 and 5, we see that the asymptotic extreme values are greater than those from the simulation of the correlated and uncorrelated series. The differences are small for the maximum values (C 1“L,) and grow to no more than 4 y0 for the seventh-highest values. The standard deviations from the asymptotic approximation are generally less than those from simulation. Compared to the simulated correlated series, the differences are not much ( 10%) for the maximum values but grow to as much as 6O”/dfor the seventh-highest values. From the foregoing, we see through simulation of variance trend

of the series

is identical

to the

deviations (in parenthexs) of the seven highest approximation of trended, uncorrciated series

Year

Max

2nd

3rd

4th

5th

6th

7th

1969

4.441 (0.182) 4.380 (0.161) 4.253 (0.174) 4.133 (0.201) 4.151

4.303 (0.110) 4.252 (0.098) 4.115 (0.104) 3.914 (0.131) c”o:E,

4.231 (0.082) 4.188 (0.072) 4.046 (0.079) 3.894 (0.104) 3.931

4.183 (0.066) 4.145 (0.057) 3.999 (0.065) 3.841 (0.090) 3.882

4.147 (0.055) 4.113 (0.048) 3.965 (0.072) 3.802 (w8:)

4.118 (0.046) 4.088 (0.040) 3.937 (0.067) 3.770 (w;)

4.094 (0.040) 4.066 (0.034) 3.914 (0.063) 3.743 w;)

3.816 (0.098) (:%)

(0.097) 3.756 (0.079) (::!:)

(0.084) 3.717 (0.069) (8%)

(0.070) 3.664 (0.057) 3.695 (0.072) 3.610 ‘;:;;;)

(0.065) 3.644 (0.054) 3.670 (0.068) 3.642 (;:;;;)

(0.069) 3.824 (0.067)

(0.064) 3.801 (0.063)

1970 1971 1972 1973

1976

(0.186) 3.934 (0.154) 4.041 (0.192) 4.047

3.882

3.799

1977

(0.210) 4.028

(0.136) 3.885

(0.108) 3.813

3.744 (0.093) 3.766

(0.076) 3.687 (0.061) 3.725 (0.078) 3.703 (0.083) 3.730

1978

(0.182) 4.14i

(0.119) 4.003

(0.176)

(0.108)

(0.095) 3.933 (0.092)

(0.082) 3.887 (0.080)

(0.074) 3.852 (0.073)

1974 1915

1721

Statistics of extreme values of air quality

the trended series that the presence of autocorrelation reduces the extreme values and increases their standard deviations. For the average extreme values, the relative differences between the correlated and uncorrelated series decrease as the rank of interest becomes less extreme; the opposite is generally true, however, for the standard deviations of the average extreme values even though the magnitudes of the standard deviations decrease by a factor of typically two to three from the maximum value to the seventh-highest value. For the extreme values and their standard directions, the asymptotic results deviate even more than the simulation results of trended uncorrelated series, when compared with the simulation results of trended correlated series. The accuracy of both the average values and their standard deviations deteriorates as the rank of interest moves away from the upper extreme. It is also of interest to compare the extreme-value estimates of different approaches using a standardized scale and without the influence of the deterministic trend. Such a comparison is presented in Appendix B. We have been dealing with the logarithms of the oxidant concentrations (in pphm). Since a small difference between two relatively large logarithmic values (typified by the extreme values) can become a significant difference in concentration space, it is worthwhile to determine what the predicted extreme concentrations are based on the asymptotic results and on the time series simulation. To determine these concentrations and their standard deviations from simulation for a specific year, we simply collected the 100 values for each rank from the simulated trended series obtained previously, took their antilogarithms and calculated the averages and standard deviations. For the asymptotic case, the means and variances of the extreme concentrations were determined by setting f(x) of (2) equal 1to eXand (eX- FqJ2, respectively and

using (3) as the cumulative distribution. This approach is straightforward. There is no need to introduce an extra approximation a/a 4 cc (Singpurwalla, 1972) just to acquire the double exponential form for G&(x). In our application, a/a is between 6 and 7.5. Table 7 shows the comparison between the asymptotic result and the simulation of the trended COTrelated series. The units are in pphm. Notice that the mean maximum value for the asymptotic case can be significantly higher than those from simulation. The overestimates get worse as the rank m increases. Except for the maximum values, the standard deviations for the asymptotic case are small compared to the simulation results. Again, the differences increase with m. For the maximum values, inclusion of the trend must have caused a greater increase in the variance of the asymptotic case than of the simulation case. The above features are qualitatively similar to the logarithmic results shown in Tables 4 and 6. Note that for the average maximum value for 1969, a 2% difference between the two estimates in the logarithmic space becomes an 11% difference in the concentration space. Compared to the single-series observed values, the simulated values are 23-72% (typically 40%) higher for the maximum concentrations, and 14-33 % (typically 20%) higher for the seventh highest concentrations. These results are comparable to those obtained using a direct lognormal fit to the original data. The consistency of overestimation clearly spells trouble for the utility of the two-parameter lognormal distribution for estimating extreme concentrations of air quality. 4. DLSCUSSlON By means of simulation we have performed a detailed study of extreme values of air quality time

Table 7. Averages and standard deviations (in parentheses) of the seven highest estimated oxidant concentrations (ii pphm) Asymptotic approximation 5th 3rd 4th

Year

Max

2nd

1969

87.0 (18.3) 81.0 (14.9) 71.5 (14.4) 63.8 (15.1) 64.8 (14.0) 51.7 (8.7) 58.1 (13.0) 58.6 (14.5) 57.2 (12.0) 64.0 (13.0)

74.5 (9.1) 70.6 (7.6) 61.6 (7.2) 53.6 (7.3) 55.2 (6.9) 45.6 (4.5) 49.3 (6.4) 49.0 (6.9) 49.0 (6.0) 55.1 (6.5)

1970 1971 1972 1973 1974 1975 1976 1977 1978

69.1 (6.4) 66.1 (5.4) 57.4 (5.2) 49.4 (5.1) 51.2 (4.9) 42.9 (3.3) 45.5 (4.5) 44.9 (4.8) 45.5 (4.2) 51.3 (4.6)

Simulation 4th

6th

7th

Max

2nd

3rd

78.3 (14.0) 19.6 (14.6) 68.2 (14.0) 60.7 (13.1) 62.6 (14.8) 50.0 (8.5) 55.1 (13.1) 53.9 (13.5) 55.5 (11.8) 58.4 (12.1)

61.9 (9.4) 68.8 (10.6) 59.7 (9.9) 51.3 (8.8) 53.3 (9.2) 44.4 (6.4) 47.8 (8.9) 46.6 (9.0) 47.9 (7.9) 50.2 (7.0)

62.0 (8.0) 63.2 (8.2) 55.0 (8.7) 46.3 (6.6) 48.4 (7.0) 41.1 (4.8) 43.4 (6.1) 42.1 (6.9) 43.1 (6.0) 46.2 (5.3)

65.8 (5.1) 63.3 (4.4) 54.7 (4.1) 46.8 (4.1) 48.7 4(;;)

63.4 (4.4) 61.3 (3.8) 52.8 (3.5) 44.9 (3.4) 46.9 2;)

61.6 (3.9) 59.7 (3.3) 51.4 (3.1) 43.5 (3.0) 45.5 J;;)

60.1 (3.5) 58.4 (3.0) 50.2 (2.8) 42.3 (2.7) 44.4 ;t;)

(2:7) 43.2 (3.6) 42.4 (3.9) 43.3 (3.4) 48.9 (3.7)

(2:3) 41.6 (3.0) 40.7 (3.2) 41.8 (2.9) 47.2 (3.2)

(2:0) 40.3 (2.7) 39.4 (2.8) 40.6 (2.5) 45.9 (2.8)

(1.8) 39.3 (2.4) 38.3 (2.6) 39.6 (2.3) 44.8 (2.5)

58.1 (6.7) 59.7 (6.5) 51.8 (7.6) 43.4 (6.1) 45.3 g.;)

.

(4.2) 41.0 (5.7) 39.0 (6.2) 41.2 (5.4) 43.8 (5.0)

5th

6th

7th

55.3 (5.7) 56.4 (6.1) 49.0 (6.5) 40.7 (5.0) 42.8 (4.2) 37.7 (3.6) 38.9 (5.3) 36.7 (5.6) 39.2 (5.2) 41.6 (4.6)

53.3 (5.6) 53.6 (5.3) 46.7 (5.7) 38.8 (4.4) 41.1 (3.8) 36.6

51.2 (4.9) 51.6 (4.8) 44.7 (5.0) 37.2 (4.0) 39.5 (3.5) 35.4

37.2 (3.5) (4.8) 35.0 37.4 (5.0)

3(Z) (4:6) 33.5 (4;)

(4.0) 40.0 (4.1)

(3.6) 38.6 (3.9)

1722

DabID P. CHOCK

series using actual data. Our approach is applicable to a very general class of time series and is devoid of the major pitfalls or statistical inconsistencies present in many popular approaches for estimating air-quality extreme values. The major assumption in our simulation is that the detrended logarithmic series fits a Box-Jenkins model with the residual series being white noise. For a series with 122 observations, the mean extreme values from simulation are relatively independent of the model used. The idiosyncrasy of the model is manifested more readily in the standard deviations. Simulation reveals that for a de&ended series with 122 observations, the asymptotic distribution gives a relatively accurate mean maximum value and standard deviation when p1 is small (say, C 0.4) but overestimation of the mean and underestimation of the standard deviation become increasingly evident as p 1 increases. The asymptotic results deteriorate esen more as the rank m increases. If the number of observations increases to 365, the performance of the asymptotic extreme-value distributions relative to simulation is expected to improve. Unfortunately, however, other undesirable and serious problems creep in as have been alluded to earfier, i.e. the problem of mixing data of different characteristics even after the seasonal trend is removed, the problem of increased ambiguity in the choice of the parental distribution, and the increased dependence of the extreme values on this distribution. When the trends were added back to the detrended series, both the extreme values and their standard deviations are perturbed. The extreme values are increased by more than the mean of the trend but less than the highest value of the trend. The standard deviations are enhanced in general, especially for the maximum values, where the enhancement is more prominent in the asymptotic case than in the simulated case. When the logarithmic results are translated back to concentration space while retaining the lognormal assumption for the uncorrelated variates in the detrended series, it is evident that si~ifi~nt relative and absolute differences exist between the extreme-value estimates from simulation and those from asymptotic distributions. The relative differences increase as the rank m increases. Even though the absolute differences are expected to decrease as the mean extreme values decrease, one must conclude that the usefulness of rhe asymptotic approxi~tion must be considered Iimited when the air quality series contains relatively high autocorrelation, as in the case with oxidant data, and when the number of observations is, say, no more than 150. The normal assumption for the white noise in the logarithmic space is largely responsible for the consistently high extreme-value estimates over those of the single-series observed data. The problem with the lognormal distribution is that it has a long tail in the upper extreme. This is compounded by the fact that a

slight absolute dit’ference in the upper extreme in the logarithmic space is magnified in the concentration space by a factor which increases with the concentration value. The problem of overesCmation becomes even worse as the extreme value increases, by increasing the number ofobservations in the series to say 355, and/or by considering very small m only. If a less extreme value is considered and if the number of 0bservatior.s is ‘kept at say. 122, the high-pollution season only, then the long-tail effect of lognormal distribution is mitigated. For example, the simulationestimates of the seventh-highest values are typically 20”; higher than the single-series observed data, compared to 400,; for the maximum values. An important additional benefit in choosing a larger m value is the reduction in the uncertainty of the estimates. For the seventh-highest value, the typical uncertainty associated with the estimate is 10 O0of the estimate, compared to 20 “0 for the maximum value. In terms of absolute concentration, the latter uncertainty is typically three times greater than the former uncertainty. From the above, we can perhaps appreciate the difficulty and the uncertainty involved in estimating the extreme values of realistic air quality data. If an air quality standard is to be tied to a concentration estimate that has the minimum uncertainty, then one ought to exclude the use of extreme values. or perhaps explore the possibility of a near-extreme standard. Of course, this standard can be made as stringent as necessary by adjusting rhe standard level. Such a standard will be less depndent on rhe serial correlation, the trend, and the choice of a fitted distribution for estimating the value of the rank of interest. In fact, even the observed value can be used directly for comparison with the standard with a relatively high level of confidence. Ackno&dgemenr-The author would like to thank Paui S. Sluchak for doing the numerical computations. REFERENCES Box G. E. P. and Jenkins G. ?A. (1976) Time Series Anai&r: ~o~ee~ring and Conrrol. Holden-Day. San Francisco,

California. Chock D. P. (1982) On the nonrandomness of high-pollution days. Atmospheric Encironment 16, 2835-2862. Chock D. P. f1984) Statistics ofextreme values of a first-order Markov normal process-an exact result. Armospheric Encironment 18. 2461-1470. Chock D. P., Terre11T. R. and Levitt S. 3. t 1975) Tie-series analysis of Riverside, California air quality data. Atmbspheric Enrironment 9.978-989. Cleveland W. S. I 19721The inverse autocorrelations of a time series and their applications. T’chnomerrics 14, 237-293. Hirtzel C. S. and Quon J. E. (1981) Statistical analysis of continuous ozone measurements. Armospheric Encironmen? 15. 1025-1034. Hirtzei C. S., Corotis R. B. and Quon J. E. l1981)Estimating the maximum value of autocorrelated air quality measurements. Atmospheric Encironmenc 16. X03-2608. Hirtzel C. S. (1985) Discussion: Statistics oiexrremc values of a first-order Markov normal process-an exact result. Atmospheric Enrironmenr 19. 1?07-1X+.

1723

Statistics of extreme values of air quality Horowitz J. (1980) Extreme values from a nonstationary stochastic process: an application to air quality analysis. Technomerrics 22. 46e478. Horowitz J. and Barakat S. (1979) Statistical analysis of the maximum concentration of an air pollutant: effects of autocorrelation and nonstationary. Armospheric Encironmenr 13,81 I-818. IMSL (1980) IMSL Library Refirence Manual (8th Edn) International Mathematical and Statistical Library, Inc., Houston, Texas. Johnson T. R. and Symons M. J. (1980) Extreme values of Weibull and lognormal distribution fitted to ambient air quality data. 73rd annual meeting of the Air Pollution Control Association, Montreal, Quebec, 22-27 June. Kendall M. (1976) Time-Series. Hafner Press, New York. Mihram G. A. (1972) Simulation: Statistid Foundncions and Melhodolo&. Academic Press, New York. Roberts E. M. (1979) Review of statistics of extreme values with applications to air quality data--?. Applications. J. Air Pollur. Control Ass. 29. 733-740. Singpurwalla N. D. (1972) Extreme values from a lognormal law with applications to air pollution problems. Technometrics 14, 703-71 I. APPESDIX

has only one peak at lag 1. This points to a (2,0.0) model with 4, positive and #t negative. The 1976 series has an inverse autocorrelation that changes signs from one lag to the next. which indicates that 4 is unlikely to be zero. The autocorrelation has two nonvanishing peaks at lags 1 and 2, the latter being much less than the square of the former. After several trials, a second-order moving-average model looks most satisfactory. The 1977 series shows one peak at lag 1 in the autocorrelation. two decaying peaks with opposite signs at lags 1 and 2 in the inverse autocorrelation. A first-order moving average model fits very well. For 1978, the autocorrelation decays essentially exponentially at small lags. Each oi the inverse and partial autocorrelations has only one peak, at lag 1. A (1.0.0)model should do well. But the autocorrelation also has a negative region at lags around 18. Inclusion of a seasonal autoregressive parameter of period I8 improves the model fit significantly. The models. together with their parameters are shown in Table 2. The chi-square values obtained from the white-noise tests on the residual autocorrelations correspond to a signitionce level of 0.4 or more. An exception is the 1973 model which gives a significance level ofO.115. The typical number of degrees oi ireedom is 21.

A

For 1969, the autocorrelation of the detrended series is prominent only at lag 1, whereas the inverse autocorrelation decays gradually with the sign alternating from one lag to the next. One immediately suspected q = 1 with 0, being negative. The autocorrelation for the 1970 series has a distinct negative peak at lag 1I. A weak negative peak also occurs at the same lag in the partial autocorrelation. Thus, a seasonal model was chosen. We do not try to rationalize why there is a period of 11 in the series. Our aim is to try to reproduce the detrended series as faithfully as possible before introducing the white noise. The 1971 series also shows a very weak periodic behavior at s = 12. The 1972 series shows a very typical first-order autoregressive behavior. The 1973 series has an autocorrelation that decays gradually with mild oscillation. It turns negative at lag 4. The inverse autocorrelation has only one spike at lag 1. The partial autocorrelation has two spikes: at lags 1 and 4. After trying several models, only the fourth-order autoregressive model was found to be satisfactory. For 1974, a (1.0.0) model is strongly indicated even though there is a mild oscillation, within + 2 standard deviations about zero in the autocorrelation. Fitting a higberorder autoregressive model showed that only 4, was significant. The 1975 series has a decaying oscillatory autocorrelation. the partial autocorrelation has significant peaks of opposite signs at lags 1 and 2, and the inverse autocorrelation

APPENDIX

B

It is of interest to compare the extreme-value estimates of different approaches using a standardized scale and without the influence of the deterministic trend. For this purpose, the series means were set equal to zero and the detrended series variances were scaled to unity. Included in the comparison, in addition to the simulation results. are the exact results based on the assumption of a first-order Markov process (Chock. 1984) and the asymptotic results for uncorrelated series. It should be noted that the process described in Chock (1984) is actually a two-state process whose transition probabilities from one state to the other is a function of autocorrelation. This process is also referred to by Hirtzel et al. (1982, see also Hirtzel. 1985) as a pseudo-Markov process. Table Bl shows a comparison of extreme-value estimates based on this two-state process assumption (only the first four ranks are available), simulation and asymptotic approximation. Our comparison is restricted to first-order autoregressive series. Let us first compare the simulation and exact results. From Table Bl. the mean extreme-values of both approaches are within 2”/, of each other for each rank. The agreements for the cases with p, > 0 are actually better than

Table Bl. Comparison ofaveragesand standard deviations (in parentheses) in standardized units ofextreme values from selected detrended series using ditferent estimation approaches

p, =

0.401

series

p, =

0.626

Exact Istorder Markov process

p, =

0.401

p, =

0.626

Simulated

( 1,w

Simulated Uncorrelated Exact Series Asymptotic

Max

2nd

3rd

4th

5th

6th

7th

2.531 (0.399) 2.468 (0.443)

2.230 (0.339) 2.211 (0.379)

2.040 (0.296) 2.010 (0.328)

1.902 (0.271) 1.896 (0.314)

I.793 (0.243) 1.788 (0.316)

1.705 (0.227) 1.701 (0.300)

1.630 (0.214) I.640 (0.298,

2.544 (0.436) 2.494 (0.445)

2.203 (0.328) 2.183 (0.369)

2.009 (0.283) 1.995 (0.324)

1.870 (0.257) 1.859 (0.297)

2.618 (0.452) 2.570 (0.422) 2.624 (0.415)

2.267 (0.322) 2.218 (0.303) 2.302 (0.260)

2.054 (0.261) 2.021 (0.255) 2.140 (0.219)

1.924 (0.246) 1.880 (0.228) 2.033 (0.172)

1.790 (0.218)

1.686 (0.196)

1.611 (0.175)

1.952 (0.152)

1.888 (0.138)

1.831 (0.126)

1723

DAVID

the case with p, = 0 where the two results should converge. The standard deviations of the simulation results are generally in good agreement with the exact results, with a tendency of being slightly greater (by no more than 8 9;) than the exact results. In terms of variances, these differences are within the uncertainties of the sample variances of the extreme values of the simulated series. It appears then that the exact results of the Xfarkov process represent the behavior of the extreme values of a first-order autoregressive series very well. In this connection, we have also compared the exact results with the results of other simulated time series of different models which have comparable p, values. While the mean extreme values are within no more than 3 Y’ (mostly less than 2 O,,)of each other for each rank, the standard deviations differ more substantially, sometimes by as much as 18 ‘& However, the general behavior of extreme values for these different time

P.

CHOCK

series remain the same. Namely, they decrease, accompanied by increasing uncertainty, as p, increases. From Table Bl. the asymptotic approximation overpredicts the average maximum value and underpredicts its standard deviation, both by a small amount, compared to the exact as well as simulated results for uncorrelated series. But the approximation deteriorates as the autocorrelation p, increases and as the rank moves away from the extreme. The mean asymptotic maximum value is 2.3 and 5 X, greater than the exact result for p, = 0 (uncorrelated series), 0.401 and 0.626, respectively, while its standard deviation is 2 ?.. 5 ?$, and 7 5; less. The mean asymptotic fourth-highest value is 8,9 and 9Y;,greater than theexact results for the threept’squoted above, and its standard deviation is 25, 33 and 42 P; less. In short, the asymptotic approximation tends to give too much confidence to the mean extreme values which it overestimates.