evaluation and comparison of statistical forecast models for daily ...
Atmospheric EnvironmentVol.24B,No. 2, pp. 303 312,1990.
0957-1272/90$3.00+0.00 © 1990PergamonPressplc
Printedin GreatBritain.
EVALUATION A N D COMPARISON OF STATISTICAL FORECAST MODELS FOR DAILY MAXIMUM OZONE CONCENTRATIONS S. M. ROaESON* a n d D. G. STEYN Department of Geography, University of British Columbia, Vancouver, British Columbia, Canada
(First received 6 December 1988 and in final form 15 June 1989)
Abstract--Three statistical models that estimate daily maximum ozone (03) concentrations in the lower Fraser Valley of British Columbia (BC) are specified using measured concentrations from two monitoring stations during the time period 1978-1985. The three models are (I) a univariate deterministic/stochastic model, (2) a univariate autoregressive integrated moving average (ARIMA) model, and (3) a bivariate temperature and persistence based regression model. The three models as well as a persistence forecast are tested by comparison with 03 concentrations observed during 1986;it is concluded that the bivariate model is superior to both univariate models and persistence. The ARIMA model has nearly the same predictive capability as persistence while the mixed deterministic/stochasticmodel performs the worst. This suggests that the traditional time series technique of decomposing a series into a trend, a cycle and a stochastic component may not be appropriate for 0 3 air quality forecasting. Key word index: Ozone, forecasting, statistical models, model evaluation.
1. INTRODUCTION
2. DATAAND STUDYREGION
Near the earth's surface, 03 is a highly toxic and reactive pollutant. Damage to human health (Goldsmith and Nadel, 1969; Bates et al., 1972; Mustafa and Lee, 1976), vegetation (Heck et al., 1983; Reich and Amundson, 1985) and materials (OECD, 1979) has been clearly demonstrated. In order to implement management and public warning strategies for 03 levels (particularly near densely populated areas), accurate forecasts of the atmospheric concentration of 0 3 are necessary. Many previous investigations have presented forecasting models for daily maximum 0 3 concentrations (e.g. McCollister and Wilson 1975; Aron and Aron, 1978; Wolffand Lioy, 1978; Prior et al., 1981; Simpson and Layton, 1983); however, very little research has been devoted to comparing existing methods to determine their relative utility. This paper uses 0 3 concentrations observed in the lower Fraser Valley of BC to develop three statistical forecast models of daily maximum 0 3 concentration. As a benchmark, persistence from one day to the next also is used as a forecast. All models are evaluated and compared employing data not used in model development.
The lower Fraser Valley region of BC, Canada, extends from the Strait of Georgia in the west to the Fraser canyon in the east. The valley is bounded by the Coast Mountains to the north and the Cascade Range to the southeast. Previous studies of atmospheric pollutants in the lower Fraser Valley by Concord Scientific Corporation (CSC, 1982, 1985) have indicated that 03 concentrations within the valley often have exceeded the maximum acceptable level (82 ppb/v) and occasionally have exceeded the maximum tolerable level (153 ppb/v) as defined by the Canadian National Air Quality Objectives. CSC (1982) determined that the largest source of the precursor pollutant emissions that lead to 03 formation in the lower Fraser Valley are motor vehicles although stationary sources (i.e. commercial and industrial fuel combustion, gasoline refining and retailing, etc.) were found to be important as well. The most densely populated and industrial area of the lower Fraser Valley lies in the Greater Vancouver Regional District (GVRD, population of approximately 1.5 million); therefore, most emissions are believed to originate there. When compared with other major urban areas in Canada, the quantity of emissions in the lower Fraser Valley is not particularly great (Dickson and Quickert, 1975); however, 0 3 concentrations reach very high levels due largely to the physiography of the area. The mountains act as a barrier to air pollutants,
*Present affiliation: Center for Climatic Research, Department of Geography, University of Delaware, Newark, DE 19716, U.S.A. 303
304
S.M. ROBESONand D. G. STEYN
channeling them along the valley. Like many coastal regions, the lower Fraser Valley experiences a seabreeze circulation that advects cool marine air, thereby limiting the depth of the planetary boundary layer (PBL) (Steyn and Oke, 1982). Together, the topography and shallow PBL combine to limit the volume of atmosphere into which pollutants may disperse. Persistent summer anticyclones intensify the 03 problem by promoting environmental conditions conducive to the production and accumulation of photochemical pollutants: warm air temperatures, light winds and strong solar input. Data used in this study are drawn from a network of air quality and meteorological monitoring stations in the lower Fraser Valley of BC operated by the GVRD and the British Columbia Ministry of Environment and Parks. A digital archive of selected pollutant and meteorological variables was established by CSC (1982) and is now available from 1978 onward; the Pollution Control Section of the GVRD updates the archive on a monthly basis. For more information regarding monitoring practices and availability of data in the lower Fraser Valley, Robeson (1987) may be consulted. In order to develop forecast models of the temporal variability of 0 3 in the lower Fraser Valley, 8 years
(1978-1985) of data from two monitoring stations were analyzed. Station T9 in Rocky Point Park and station T l l at Abbotsford Airport (see Fig. 1) were chosen on the basis of the range of historical 03 concentrations, their spatial location and separation, and the quality and length of data records. The two stations are aligned roughly parallel to the wind directions (i.e. westerly sector) that commonly produce the highest 0 3 concentrations (CSC, 1985). In addition, they are crudely representative of two distinct atmospheric situations--an urban or 'local transport' situation (station T9) and a rural or 'regional transport' situation (station T11). Local transport is characterized by strong local emissions causing high pollutant concentrations near emission sources. In contrast, regional transport involves emissions that travel far from their sources and may undergo extensive chemical transformation before reaching a downwind location. Station T9, located in Rocky Point Park alongside Burrard Inlet in the Port Moody Basin, was established in 1977 to monitor general community air quality in the city of Port Moody. The Environmental Protection Service of Environment Canada has designated T9 as a National Air Pollution Survey Class 1 station and contributes financial and technical sup-
Fig. 1. Map of the lower Fraser Valley, BC study area showing the topography, coastline, built-up areas and the two monitoring stations used in the study.
Statistical forecast models for ozone port. Just to the west of the station are wood processing and coal/sulphur loading terminals while farther to the west are four petroleum refineries and an intermittently-used natural gas-fired power generating plant. Given the large local sources of emissions, the heavily populated areas to the south and west, and the potential for local stagnation of airflow, there is obviously great concern regarding air quality in the Port Moody Basin. Ozone concentrations at station T l l (Abbotsford Airport) have also been found to exceed air quality standards (CSC, 1985). In contrast to station T9, station T11 is located near a lightly-used airport in an area of largely agricultural land-use where no large local sources of emissions exist. Station T 11 is directly downwind of the major emissions source area (i.e. the greater Vancouver region) when O3 concentrations throughout the valley are highest (CSC, 1985). Hence, it is believed that the source of the photochemical pollutants that reach station T11 is the greater Vancouver region. 3. D E S C R I P T I O N
305
In order to reproduce the seasonal variation of Oa concentrations, Horowitz and Barakat fitted a second-order polynomial (regressed on Julian day number t and on t 2) to data that had first been logtransformed. The use of a polynomial to describe annual variability appears appropriate in as much as the ordinary least squares procedure objectively determines the timing of the summer maximum. Horowitz and Barakat (1979) applied a simple log transformation to the data before determining the least squares coefficients of the polynomial; here, the data are transformed via the three-parameter lognormal procedure outline by Ott and Mage (1976). A second-order polynomial is sufficient to remove any trend in the data from stations T9 and T l l . After removal of the polynomial, a highly autocorrelated sequence of data (e,) remains: et =
l n ( x , - K ) - (b o + bit
Three forecast models are described and developed using data from stations T9 and T l l during the months of May to September, 1978-1985. Only high 03 months are chosen to avoid mixing data that have different characteristics (Chock, 1986). The comparative ability of the models to forecast daily maximum 1-h average 03 concentrations is presented in section 4.
et = dplet- 1 + at
(1)
(2)
where ~b1 is a first order autoregressive coefficient, e,_ is the value of the residual series on the previous day, and a t is a normally and independently distributed sequence with a mean value of approximately zero. Hence, all trend and serial correlation have been removed and the forecast O a concentration, ~t, is described by
3.1. A deterministic/stochastic (D/S) model Horowitz and Barakat (1979) have identified the general form of a non-stationary stochastic process that is able to simulate observed pollutant concentrations without assuming independence or first-order stationarity. They examined one year of 03 data from the St. Louis Regional Air Pollution Study to demonstrate the usefulness of such a model.
~t=exp[(bo+blt+b2t2)+dPlet_l]+K.
(3)
Parameter estimates for stations T9 and T l l using data from 1978 to 1985 are presented in Table 1 (1986 data were reserved to evaluate the model's performance). Previously, this model has been used to
Table 1. Parameter estimates (1) Deterministic/stochastic model Station bo b1 T9 3.57 8.58 x 10-2 (0.02)* (2 x 10-4) Tll 3.81 4.45 x 10-2 (0.01) (2 x 10 -4)
K
b2
-2.45 × 10 -4 (6 x 10- ~) - 1.33 x 10 -4 (5 x 10- 7)
b2
0.309 (0.04) 0.433 (0.04)
*Numbers in parentheses are standard errors of estimates.
AE(B) 24:2-[
b2 t2)
where Zt is the ozone concentration on day number t, K is a constant, and b o, b 1 a n d b 2 are ordinary least squares regression coefficients. Inspection of autocorrelation and partial autocorrelation functions for the residual series et clearly suggests a first-order autoregressive equation (for all years):
OF MODELS
(2) ARIMA (1, 1, 0) Model Station ~bl T9 0.739 (0.08)* T11 0.756 (0.08) (3) TEMPER Model Station bo b1 T9 - 16.3 2.48 (3.7)* (0.19) Tll 1.51 1.15 (2.8) (0.14)
+
0.526 (0.02) 0.592 (O.O2)
-25 -25
306
S.M. ROBESO~qand D. G. STEV~q
forecast probabilistic aspects of Oa time series (Robeson and Steyn, 1989). The model is termed deterministic/stochastic (D/S) in reference to the seasonal polynomial (deterministic component) and the autoregressive term (stochastic component).
Mean values of estimated model parameters for stations T9 and T l l are given in Table 1. The ability of the ARIMA (1,1,0) model to forecast daily maxima of O a is discussed in section 4. 3.3. A temperature and persistence based model
3.2. An A R I M A model The time series method of Box and Jenkins (1976) is a commonly used technique for forecasting O3 concentrations (e.g. Merz et al., 1972; Sandberg et al., 1974; Chock et al., 1975; McCollister and Wilson, 1975; Simpson and Layton, 1983). A purely stochastic autoregressive integrated moving average (ARIMA) model may be developed. When given an original series L, the general form of the model is represented by w , = 471 w t - 1 + "
" " + dpp w t _ t, + a t - 0 1 a t - 1
.....
Oqar_ q
(4) where w, is the differenced original series, w, = VZ,.
(5)
Box and Jenkins (1976) define the backward difference operator as V~(,= X,- Z,- 1
(6)
which may be successively applied to the series d times. A troublesome property of the differencing filter is its 'noise amplification' (Bennett, 1974). Although in practice d rarely exceeds 2 (Box and Jenkins, 1976), differencing should be used carefully since the magnification of high frequency components may be consequential. The autoregressive (qb) and moving average (0) model parameters are determined either by maximum likelihood or least squares estimation. The effects of all factors other than past time series values are incorporated into the 'random shock' term, at. Proper model identification should result in an a, series that is normally and independently distributed with a mean value of approximately zero. Inspection of autocorrelation and partial autocorrelation functions for the 0 3 series suggests that the most appropriate model (i.e. the one with relatively small white noise variance and a small number of parameters) is an ARIMA (1,1,0) model, where the notation ARIMA (p,d,q) refers to the number of autoregressive parameters, p, moving average parameters, q, as well as the number of differencing operations performed, d. Various versions of Equation (4) were tested to find a model that accurately describes seasonal variability and maintains some degree of parsimony; however, with increasingly complex models, very small decreases in residual variance resulted. An ARIMA (1,1,0) model has the form w t + q91 w , _ x = at.
(7)
When transformed to the original series, this becomes ~t = (~b 1 + 1 ) Z , _ I - 4~1Z,- 2 + at.
(8)
The two models previously described are variations of the 'pure time series' approach, i.e. a univariate series is decomposed into various elements. In this section, a model that relates daily maximum 0 3 to 'external' variables as well as to previous O3 values is developed. Forecast models often rely heavily on regression techniques for both model formulation and estimation (Aron and Aron, 1978; Wolff and Lioy, 1978; Aron, 1980; Prior et al., 1981). Numerous independent variables commonly are entered into a stepwise regression to identify those having 'significant' partial correlations with daily maximum 0 3. Such models typically explain a non-trivial proportion of the variance in daily O a maxima; however, practical and methodological difficulties preclude their use here. First, many of the independent variables commonly used are derived from upper air data (e.g. Aron and Aron, 1978). Such data were not available for this study since no upper air stations are located in the lower Fraser Valley. Second, simply entering a large number of independent variables into a stepwise regression procedure is neither physically purposeful not intuitively appealing. When selecting independent variables to forecast 0 3 variability, factors that influence both photochemical production and atmospheric accumulation should be considered. Rudimentary examination of Oa photochemistry indicates the importance of NO, NO 2 and certain species of HCs (Fishman and Crutzen, 1978). Ozone accumulation is related to dispersive properties of the atmosphere and hence to wind speed, inversion height and stability. Solar radiation (in particular, u.v. wavelengths) also affects O3 variability in a well-known fashion. Munn (1975), however, has shown that the strength of solar radiation input across Canada is not a limiting factor in photochemical reactions during summer. Of all meteorological variables, air temperature has the strongest correlation with 03 concentrations for two reasons: (1) high air temperatures are an excellent indication of environmental conditions conducive to 03 production and accumulation (i.e. anticyclonic conditions with associated clear skies and light winds) and (2) the rate 'constants' of photochemical reactions are highly temperature dependent (OECD, 1979). As a result (at least for the lower Fraser Valley), air temperature is a reasonable surrogate for the combined effects of wind speed, wind direction, inversion height and photochemical reaction rate. Consequently, correlation and partial correlation analyses were performed to identify important independent variables from a set including NO, NO2, wind speed and its inverse, u.v. radiation,
Statistical forecast models for ozone air temperature and the previous day's O 3 concentration. Speciated HC data were not available. A combined linear relationship between daily maximum 03 concentration and (I) daily maximum air temperature (Tt) and (2) the previous day's daily maximum 03 concentration (Xt-l) explains a large proportion of the variance in Xt. This may be expressed as
307
regression of forecast variable on observed variable, mean absolute and root mean squared errors (MAE, RMSE), and systematic and unsystematic components of RMSE: [-1 n
^
271/2
RMSE~=|- ~ (P,-O,) / Lni=l
(10)
J
[-1 ~ 2]1/2 RMSEu = Lni~I(PI- Pi) ] , .
f(, = bo + bl Tt + b2z,- ~
(9)
where bo, bl and b2 are coefficients from an ordinary least squares regression (see Table 1). Explained variance (r 2) at station T9 is 0.52 while for station T l l , r 2 = 0.41. These values of explained variance are slightly lower than those obtained by Prior et al. (1981) via a similar approach; however, unlike that study, the data used here (1) contained multi-year variability, (2) did not include early morning values, and (3) were not spatially averaged. As a result, forecasts based on Equation (9) are both timely and meaningful for operational purposes. Daily maximum air temperature data were taken from 1-h average air temperature values at station T9 from 1982 to 1985 since temperature data were not available at station T9 prior to 1982. If used operationally in a forecasting mode, Equation (9) requires accurate forecasts of daily maximum temperature. The only temperature forecast verification performed in the lower Fraser Valley occurs at Vancouver International Airport where forecasts for the months of June-August of 1986 had a mean absolute error of 1.27°C (Hammond, pers. comm., 1987). Further inland (i.e. away from the moderating influence of the Strait of Georgia), temperature deviations are likely to be larger; nonetheless, daily maximum air temperature forecast errors of I°C result in 03 forecast errors of only 2.5 ppb at station T9 and 1.1 ppb at station Tll--values small enough to warrant using forecast air temperature. This bivariate regression model, Equation (9), is termed TEMPER due to its reliance of Temperature and Persistence. Its ability to forecast daily maxima of 03 is examined in the following section. 4. M O D E L PERFORMANCE
Three statistical models of the temporal variability of 0 3 were described in the previous section; however, no attempt was made to compare the relative ability of each model to forecast daily maximum O a concentrations. In this section, the forecast performance of each model is evaluated via graphical and statistical comparison with data from the year 1986--data not used in model development. In addition, all models are compared with a persistence forecast (i.e. f~t= ;(t- ~), a benchmark comparison in nearly all of the atmospheric and environmental sciences. As recommended by Willmott (1981) and Willmott et al. (1985), an array of model evaluation statistics are presented: observed and forecast means and standard deviations, intercept (a) and slope (b) of a least squares
(11)
where Oi and Pi are the observed and predicted values, respectively and Pi = a+ boy The simple linear relationship RMSE 2 = R M S E 2 + R M S E 2 indicates the very useful decomposition of the total error into systematic and unsystematic elements. Also reported is the Index of Agreement (Willmott, 1981), d, defined as d= 1
X~=~(Pi- O~)~ , , , 2; ~ : 1(lPil + IO,I)
(12)
where P [ = P i - O and O ' i = O i - O , with 0 being the observed mean value. Being dimensionless and having the limits of 0.0 (indicating no agreement) and 1.0 (indicating perfect agreement), d may be viewed as a standardized (by the variability in the predictions and observations about the observed mean) measure of the mean squared error. To forecast daily maximum O a concentrations with the D/S and ARIMA models, the random term, at, is simply set equal to its expected value: zero. Forecasts using the TEMPER model employed actual, rather than forecast daily maximum air temperature values. To test the sensitivity of TEMPER to air temperature forecast errors, a series of normally distributed random numbers (with mean of 0.0 and standard deviation of 2.0) is added to observed temperatures. Air temperatures containing these random errors are then used in TEMPER to forecast daily maximum 03 concentrations. Scatterplots of observed vs forecast daily maximum 03 concentrations are presented in Fig. 2 (a)-(d) for station T9 and Fig. 3 (a)-(d) for station T11. Scatterplots for the TEMPER with error forecasts are not included since they are nearly indistinguishable from those of the TEMPER forecasts: Model evaluation statistics are summarized in Table 2 for comparative purposes. An examination of Table 2 reveals a few noteworthy facts. First, the persistence forecast was better than that of the D/S model and nearly the same as that of the ARIMA model at both stations T9 and Tll. Second, for all models, the slope and intercept of the regression line indicate that low concentrations are being overpredicted and high concentrations are being underpredicted; apparently, none of the models are sufficiently dynamic to capture the rapid variability in the 03 time series. Although nearly all statistical models underestimate variance (since they are usually some form of non-magnifying filter through data), the D/S model performs particularly poorly in this regard. Third, the random noise added to the daily maximum
308
S.M. ROBESONand D. G. STEYN (o) o/s
(b)
140.0
140.0.
O. 120.0 O. v
Q. 120.0,
o
O.
C 100.00
(~ 100.0.
80.0-
o
C
o
.o•
oo
,
C 0
60.0-
o*" o o o
.
o
o
o
,o0 o
,
o
o
O .~L C 0 ~)
o
o .o
o
o
20.0' o
o.0 0.0 " 20:0
40:0 "60:0
"80:0
0.0
" I0().0" 12().0" 140.0
20.0
0.0
Observed
Concentrotion
(ppb)
(d)
140.0-
o
80.0-
o
o
,"
d)oO o Y
"°°°
o o 6, %0 . ,
60.0-
0
C)
120.0
140,0
(ppb)
o
o
oo
o
100.0.
o o o
O C 0 U )
o
°°° o
o
•o.o
o oo
o
ooo
i
20.0
-
,
40.0
Observed
-
,
60.0
•
1
-
80.0
Concentrotion
,
100.0
-
,
120.0
-
140.0
(ppb)
o
oo
°o°
20.0. o
oo
~°° °o i
o
o o o
~
o
O U..
o.6
60.0.
O
o
. o
80.0
0 U
0,0
100.0
CL 120.0 CL
o /
c- 100.00
0.0
80.0
Persistence
z"
Q. 120.0O.
20.0-
60.0
Concentrotion
140.0 ,
z"
~ 0 b-
40.0
Observed
(c) T E M P E R
40.0 -
-'" ..o" ° °
40.0"
O LL
"~ 0
o
o o o ° o oo o
O U
0
r-
. e
o o
;o
60.0-
,
o
o
60.0.
O
Z:oO
"~)'oo
20.o O o ~ .
.-.-
b_
C
ARIMA
°° o °°°
o
0.0 0.0
" 2010
" 4010
Observed
"60:0
"8010
Concentration
1 0 0 . 0 " 12().0
140.0
(ppb)
Fig. 2. Scatterplots of forecast vs observed daily maximum ozone concentrations for the summer of 1986 at station T9 using: (a) D/S, (b) ARIMA (1,1,0), (c) TEMPER and (d) persistence models. Solid line indicates perfect prediction while the dashed line is an ordinary least squares regression estimate.
temperature data has had little effect on the forecast performance of the TEMPER model. A traditional time series approach (decomposition of the variability into trend, cycle and stochastic elements) was taken to develop the D/S model. However, the D/S model did not perform well--the naive persistence forecast is clearly better with the data used here. The D/S model contains the most linear systematic error as indicated by its RMSEs; the ratio MSEJMSE ~ 0.70 indicates that 70% of the total error is systematic. It is likely that appropriate modifications (e.g. relating the 'whitened' series, a,, to a 'whitened' air temperature series) would improve its forecasting ability. However, by explicitly including a polynomial representation of the seasonal cycle of 03,
high and low values in the 03 series are modulated. This is accomplished when the polynomial is subtracted from the current 0 3 concentration and the residual (rather than the entire value as in the ARIMA, TEMPER, and persistence forecasts) is used as the dynamic component in the forecast. These results also indicate that the utility of the D/S model as a conditional probability density function (Robeson and Steyn, 1989) may be compromized. In other disciplines such as economics where time series analysis also plays an important role, the traditional approach has been found to be inadequate in many situations (Granger, 1980). Hence, when using univariate methods, alternative strategies such as adaptive parameter estimation may be necessary to empirically model
Statistical forecast models for ozone (o) o/s 100.0 O. (3.
80.0
C 0
cO
60.0-
0
(..)
o o
~o
.-"
omo o
o
C 0
o
o
Qo
o o
U tO
o
0
°. . ,
0 b_ 0.0 0.0
0 (.)
2o.o
40.0
6o0
ao.o
0.0 0.0
~oo.o
0
20.0
(d)
TEMPER
o
"'"
oO ooo
~.o'o
40.0
eo.o
Bo.o
100.o
Persistence
100.0
I
°° ~
80.0.
80.0-
°
oo
o o 60.0-
.~
(9
60.0
g
o
0
40.0
4O.0.
20.0
20.0
.'"
o
. o.,
o
o
o o
o °
o ~ o ~ o . "1~,~i~ o°
~
o o
o o oo o , °° ~ , ~ o ° o
"~' J ' ~ ' C
O
o
oO
0
O O
oo
Observed Concentration (ppb)
C E
o° o o
.-
2O,0
C 0 "~
, ." o . -
0 L
100.0
O. CL
o
Oo o o o L ~ , / ~c oo oo o - oOO -" oo~, ~ o -
ObseFved ConcentFation {ppb)
(c)
o
o
40,0
o
20.0-
oo
60,0
o
C)
o . -~
o o
40.0- -o o
o
80,0
o
o
o
O.
o
o .-'"
ARIMA
100.0-
-'°"
o o ®o
o o
oo o
o
C
(J c
o
(b)
309
To
o -o
o
o
0
0.0 0.0
20:0 Observed
' 40.0
" 60.0 '
Concentrotion
" 80.0 '
100.0
(ppb)
0.0 0.0
2olo 4 o l o Observed
~o:o
Concentration
6olo
~oo.o
(ppb)
Fig. 3. Scatterplots of forecast vs observed daily maximum ozone concentrations for the summer of 1986 at station T11 using: (a) D/S, (b) ARIMA (1, 1,0), (c) TEMPER and (d) persistence models. Solid line indicates perfect prediction while the dashed line is an ordinary least squares regression estimate.
systems (such as the photochemical one) that vary on short timescales. To gain an understanding of why the ARIMA and persistence models performed similarly, it is useful to view the two models as filters: the original time series is transformed into an uncorrelated white noise sequence. The persistence model transforms a time series in exactly the same way as the first-order differencing filter described in section 3.2 does, while the ARIMA (1,1,0) model is simply an optimally chosen (in this case by maximum likelihood estimation) weighted average of previous values in the Oa series. Comparison of the amplitude spectra of the ARIMA(1, 1,0) and persistence filters (Fig. 4) reveals the reason for the similar forecasting ability--their amplitudes in the frequency domain are nearly identical (as one might
also suspect based on the magnitude of the autoregressive coefficients). Phase spectra of the two filters are similar as well. The widespread use of the Box-Jenkins (ARIMA) method in the atmospheric and environmental sciences (McLeod et al., 1977; Bennett, 1979) results from its ability to forecast state variables while maintaining model simplicity (the models presented here have only one parameter, ~bt, to estimate). However, the performance of such models should always be compared--from the perspective of model evaluation, not simply in terms of the residual variance in the fit--with that of persistence. Regardless of simulated forecast errors in daily maximum air temperatures, the TEMPER model appears to be superior to the other models as indicated by the array of model evaluation statistics (in
S. M. ROBESON and D. G. STEYN
310
Table 2. Model evaluation statistics Station T9
(1)
Observed mean Observed standard deviation Predicted mean Predicted standard deviation Regression line intercept slope Mean absolute error Root mean square error systematic component unsystematic component Index of agreement
Station T l l Observed mean Observed standard deviation Predicted mean Predicted standard deviation Regression line intercept slope Mean absolute error Root mean square error systematic component unsystemetic component Index of agreement
Model* (3)
(2) - -
41.3
(4)
(5)
48.1 18.3
41.9 21.0
21.9 0.632 12.9 16.2 10.2 12.6 0.821
18.7 0.560 13.2 19.7 9.3 17.4 0.743
- -
--21.1-49.0 17.5
42.7 11.4
41.7 19.3
29.7 0.313 13.2 17.3 14.5 9.3 0.678
21.0 0.500 13.6 19.3 10.5 16.2 0.733
(1)
(2)
Model (3)
(4)
(5)
39.6 9.4
41.5 15.2
41.3---16.6~ 43.6 11.3
43.0 11.6
41.5 16.5
26.0 0.328 9.8 13.6 11.3 7.6 0.686
21.5 0.484 10.7 15.5 8.5 12.9 0.722
23.3 0.476 9.7 12.2 8.8 8.5 0.792
17.9 0.570 10.4 15.3 7,1 13.5 0.752
22.7 0.636 12.3 15.6 10.8 11.2 0.832
23.8 0.480 9.3 12.0 8.9 8.1 0.798
*(1) D/S, (2) ARIMA, (3) TEMPER, (4) TEMPER with error, (5) persistence.
2.0
PURE PERSISTENCE
1.5 t-,.J
1.0
0.5
0.0 0'
;o
40 100 120 140 ' 2o 8'0 ' ' ' LINEAR FREQUENCY (CYCLES PER YEAR)
160 '
180 '
Fig. 4. Amplitude spectra of the filtering process represented by pure persistence and the ARIMA (1,1,0) models at stations T9 and Tll. particular, note the low M A E a n d high Index of Agreement). P e r f o r m a n c e (and the original fit) is s o m e w h a t p o o r e r at station T11 a l t h o u g h this simply m a y be a result of the t e m p e r a t u r e d a t a being d r a w n
from station T9. T E M P E R also forecasts reasonably well on days with high 0 3 concentrations, a n important prerequisite for a n operational forecast model. Since T E M P E R is a regression model, it requires
Statistical forecast models for ozone (1) accurate forecasts of predictor variables (i.e. daily maximum air temperature) and (2) periodic reestimation of model parameters as new data become available. All of the models appear to have some difficulty forecasting high 0 3 values; however, it should be noted that none of these models were developed to explicitly handle extreme values. Extreme value prediction is often handled in some ad hoc fashion (Prior et al., 198 I)---no such corrections were attempted here. Given the hazards of exceeding epidemiologically determined standards, though, the models' performance with high concentrations suggests that a different form of model (e.g. quadratic or multiplicative rather than a linear dependence) is worth exploring. 5. CONCLUSIONS Three statistical models and persistence have been used to forecast daily maximum 0 3 concentrations in the lower Fraser Valley of British Columbia. Through graphical and statistical comparison with data not used in model development, the relative utility of each model has been demonstrated. Of the three statistical models, a regression model (TEMPER) based on daily maximum air temperature and the previous day's 0 3 concentration provides the most accurate forecast. Incorporating physically meaningful variables such as air temperature clearly gives a regression model an advantage over univariate methods. An A R I M A model forecasts no better than pure persistence while a traditional time series decomposition (the D/S model) performs the worst with these data. Hence, given reasonably accurate forecasts for daily maximum air temperature, it is suggested that T E M PER be used for forecasts of daily maximum 0 3 concentrations in the lower Fraser Valley. The best 'rule of t h u m b ' - - a n d certainly the simplest--forecast is persistence. The models do appear to have some difficulty forecasting very high concentrations. Although ad hoc corrections may improve local forecasts, such models will likely be site specific. It is suggested that different forms of models (e.g. including quadratic or multiplicative dependencies) be developed and evaluated for general applicability. Acknowled#ements--The Pollution Control Section of the Greater Vancouver Regional District kindly provided the data used in this study. Thoughtful comments from Cort Willmott (University of Delaware) and two anonymous reviewers are appreciated. This research was supported by a grant from the Natural Science and Engineering Research Council of Canada and by a University of British Columbia Graduate Fellowship to SMR. REFERENCES
Aron R. (1980) Forecasting high level oxidant concentrations in the Los Angeles basin. J. Air Pollut. Control Ass. 30, 1227-1228. Aron R. and Aron I.-M. (1978) Statistical forecasting models:
311
II. Oxidant concentrations in the Los Angeles basin. J. Air Pollut. Control Ass. 28, 684-688. Bates D. V., Bell G. M., Burnham C. D., Hazucha M, Mantha J., Pengeily L. D. and Silverman F. (1972) Shortterm effects of ozone on the lung. J. appl. Physiol. 32, 176-181. Bennett R. J. (1974) Process identification for time series modelling in urban and regional planning. Regional Stud. 8, 157-174. Bennett R. J. (1979) Spatial Time Series. Pion, London. Box G. E. P. and Jenkins G. M. (1976) Time Series Analysis, Forecasting and Control. Holden-Day, San Francisco. Chock D. P. (1986) Statistics of extreme values of air quality--a simulation study. Atmospheric Environment 10, 1713-1724. Chock D. P., Terrell T. R. and Levitt S. B. (1975) Time-series analysis of Riverside, California air quality data. Atmospheric Environment 9, 978-989. Concord Scientific Corporation (CSC) (1982) Vancouver Oxidant Study, Air Quality Analysis 1978-1981, Prepared for Environment Canada, Environmental Protection Service. Concord Scientific Corporation (CSC) (1985) Vancouver Oxidant Study, Air Quality Analysis Update 1982-1984. Prepared for Environment Canada, Environmental Protection Service. Dickson D. R. and Quickert N. (1975) The chemical composition of photochemical air pollution. In Photochemical Air Pollution: Formation, Transport and Effects. National Research Council of Council, Publication No. NRCC 14096 of the Environmental Secretariat. Fishman J. and Crutzen P. J. (1978) The origin of ozone in the troposphere. Nature 274, 855-858. Goldsmith J. R. and Nadel J. A. (1969) Experimental exposure of human subjects to ozone. J. Air Pollut. Control Ass. 19, 329-330. Granger C. W. J. (1980) Forecasting in Business and Economics. Academic Press, New York. Heck W. W., Adams R. M., Cure W. W., Heagle A. S., Heggestad H. E., Kohut R. J., Kress L. W., Rawlings J. O. and Taylor O. C. (1983) A reassessment of crop loss from ozone. Envir. Sci. Technol. 572A-581A. Horowitz J. and Barakat S. (1979) Statistical analysis of the maximum concentration of an air pollutant: effects of autocorrelation and non-stationarity. Atmospheric Environment 13, 811-818. McCollister G. M. and Wilson K. R. (1975) Linear stochastic models for forecasting daily maxima and hourly concentrations of air pollutants. Atmospheric Environment 9, 417-423. McLeod A. I., Hipel K. W. and Lennox W. C. (1977) Advances in Box-Jenkins modelling 2. Applications. Water Resour. Res. 13, 577-586. Merz P. H., Painter L. J. and Ryason P. R. (1972) Aerometric data analysis--time series analysis and forecast and an atmospheric smog diagram. Atmospheric Environment 6, 319-342. Munn R. E. (1975) The oxidant climatology of Canada. In Photochemical Air Pollution: Formation, Transport, and Effects. National Research Council of Canada, Publ. No. NRCC 14096 of the Environmental Secretariat. Mustafa M. G. and Lee S. D. (1976) Pulmonary biochemical alterations resulting from ozone exposure. Ann. Occup. Hyg. 19, 17-26. Organization for Economic Co-operation and Development (OECD) (1979) Photochemical Oxidants and Their Precursors in the Atmosphere. OECD, Paris. Ott W. R. and Mage D. T. (1976) A general purpose univariate probability model for environmental data analysis. Comp. Oper. Res. 3, 209-216. Prior E. J., Schiess J. R. and McDougal D. S. (1981) Approach to forecasting daily maximum ozone levels in St. Louis. Envir. Sci. Technol. 15, 430-436.
312
S.M. ROBESON and D. G. STEYN
Reich P. B. and Amundson R. G. (1985) Ambient levels of ozone reduce net photosynthesis in tree and crop species. Science 230, 566-570. Robeson S. M. (1987) Time series analysis of surface layer ozone in the lower Fraser Valley of British Columbia. M.Sc. Thesis, The University of British Columbia, Vancouver, British Columbia, Canada. Robeson S. M. and Steyn D. G. (1989) A conditional probability density function for forecasting ozone air quality data. Atmospheric Environment 23, 689-692. Sandberg J. S., Robinson L. H. and DeMandel R. E. (1974) Box-Jenkins stochastic modeling in regional forecasting of photochemical oxidant. In AMS Conference on Atmospheric Turbulence and Diffusion, 1974. American Meteorological Society, Boston.
Simpson R. W. and Layton A. P. (1983) Forecasting peak ozone levels. Atmospheric Environment 17, 1649-1654. Steyn D. G. and Oke T. R. (1982) The depth of the daytime mixed layer two coastal sites: a model and its validation. Boundary-Layer Met. 24, 161-180. Willmott C. J. (1981) On the validation of models. Phys. Geog. 2, 184-194. Willmott C. J., Ackleson S. G., Davis R. E., Feddema J. J., Klink K. M., Legates D. R., O'Donnell J. and Rowe C. M. (1985) Statistics for the evaluation and comparison of models. J. geophys. Res., 90, 8995-9005. Wolff G. T. and Lioy P. J. (1978) An empirical model for forecasting maximum daily ozone levels in the northeastern U.S.J. Air Pollut. Control Ass. 28, 1034-1038.
0957-1272/90$3.00+0.00 © 1990PergamonPressplc
Printedin GreatBritain.
EVALUATION A N D COMPARISON OF STATISTICAL FORECAST MODELS FOR DAILY MAXIMUM OZONE CONCENTRATIONS S. M. ROaESON* a n d D. G. STEYN Department of Geography, University of British Columbia, Vancouver, British Columbia, Canada
(First received 6 December 1988 and in final form 15 June 1989)
Abstract--Three statistical models that estimate daily maximum ozone (03) concentrations in the lower Fraser Valley of British Columbia (BC) are specified using measured concentrations from two monitoring stations during the time period 1978-1985. The three models are (I) a univariate deterministic/stochastic model, (2) a univariate autoregressive integrated moving average (ARIMA) model, and (3) a bivariate temperature and persistence based regression model. The three models as well as a persistence forecast are tested by comparison with 03 concentrations observed during 1986;it is concluded that the bivariate model is superior to both univariate models and persistence. The ARIMA model has nearly the same predictive capability as persistence while the mixed deterministic/stochasticmodel performs the worst. This suggests that the traditional time series technique of decomposing a series into a trend, a cycle and a stochastic component may not be appropriate for 0 3 air quality forecasting. Key word index: Ozone, forecasting, statistical models, model evaluation.
1. INTRODUCTION
2. DATAAND STUDYREGION
Near the earth's surface, 03 is a highly toxic and reactive pollutant. Damage to human health (Goldsmith and Nadel, 1969; Bates et al., 1972; Mustafa and Lee, 1976), vegetation (Heck et al., 1983; Reich and Amundson, 1985) and materials (OECD, 1979) has been clearly demonstrated. In order to implement management and public warning strategies for 03 levels (particularly near densely populated areas), accurate forecasts of the atmospheric concentration of 0 3 are necessary. Many previous investigations have presented forecasting models for daily maximum 0 3 concentrations (e.g. McCollister and Wilson 1975; Aron and Aron, 1978; Wolffand Lioy, 1978; Prior et al., 1981; Simpson and Layton, 1983); however, very little research has been devoted to comparing existing methods to determine their relative utility. This paper uses 0 3 concentrations observed in the lower Fraser Valley of BC to develop three statistical forecast models of daily maximum 0 3 concentration. As a benchmark, persistence from one day to the next also is used as a forecast. All models are evaluated and compared employing data not used in model development.
The lower Fraser Valley region of BC, Canada, extends from the Strait of Georgia in the west to the Fraser canyon in the east. The valley is bounded by the Coast Mountains to the north and the Cascade Range to the southeast. Previous studies of atmospheric pollutants in the lower Fraser Valley by Concord Scientific Corporation (CSC, 1982, 1985) have indicated that 03 concentrations within the valley often have exceeded the maximum acceptable level (82 ppb/v) and occasionally have exceeded the maximum tolerable level (153 ppb/v) as defined by the Canadian National Air Quality Objectives. CSC (1982) determined that the largest source of the precursor pollutant emissions that lead to 03 formation in the lower Fraser Valley are motor vehicles although stationary sources (i.e. commercial and industrial fuel combustion, gasoline refining and retailing, etc.) were found to be important as well. The most densely populated and industrial area of the lower Fraser Valley lies in the Greater Vancouver Regional District (GVRD, population of approximately 1.5 million); therefore, most emissions are believed to originate there. When compared with other major urban areas in Canada, the quantity of emissions in the lower Fraser Valley is not particularly great (Dickson and Quickert, 1975); however, 0 3 concentrations reach very high levels due largely to the physiography of the area. The mountains act as a barrier to air pollutants,
*Present affiliation: Center for Climatic Research, Department of Geography, University of Delaware, Newark, DE 19716, U.S.A. 303
304
S.M. ROBESONand D. G. STEYN
channeling them along the valley. Like many coastal regions, the lower Fraser Valley experiences a seabreeze circulation that advects cool marine air, thereby limiting the depth of the planetary boundary layer (PBL) (Steyn and Oke, 1982). Together, the topography and shallow PBL combine to limit the volume of atmosphere into which pollutants may disperse. Persistent summer anticyclones intensify the 03 problem by promoting environmental conditions conducive to the production and accumulation of photochemical pollutants: warm air temperatures, light winds and strong solar input. Data used in this study are drawn from a network of air quality and meteorological monitoring stations in the lower Fraser Valley of BC operated by the GVRD and the British Columbia Ministry of Environment and Parks. A digital archive of selected pollutant and meteorological variables was established by CSC (1982) and is now available from 1978 onward; the Pollution Control Section of the GVRD updates the archive on a monthly basis. For more information regarding monitoring practices and availability of data in the lower Fraser Valley, Robeson (1987) may be consulted. In order to develop forecast models of the temporal variability of 0 3 in the lower Fraser Valley, 8 years
(1978-1985) of data from two monitoring stations were analyzed. Station T9 in Rocky Point Park and station T l l at Abbotsford Airport (see Fig. 1) were chosen on the basis of the range of historical 03 concentrations, their spatial location and separation, and the quality and length of data records. The two stations are aligned roughly parallel to the wind directions (i.e. westerly sector) that commonly produce the highest 0 3 concentrations (CSC, 1985). In addition, they are crudely representative of two distinct atmospheric situations--an urban or 'local transport' situation (station T9) and a rural or 'regional transport' situation (station T11). Local transport is characterized by strong local emissions causing high pollutant concentrations near emission sources. In contrast, regional transport involves emissions that travel far from their sources and may undergo extensive chemical transformation before reaching a downwind location. Station T9, located in Rocky Point Park alongside Burrard Inlet in the Port Moody Basin, was established in 1977 to monitor general community air quality in the city of Port Moody. The Environmental Protection Service of Environment Canada has designated T9 as a National Air Pollution Survey Class 1 station and contributes financial and technical sup-
Fig. 1. Map of the lower Fraser Valley, BC study area showing the topography, coastline, built-up areas and the two monitoring stations used in the study.
Statistical forecast models for ozone port. Just to the west of the station are wood processing and coal/sulphur loading terminals while farther to the west are four petroleum refineries and an intermittently-used natural gas-fired power generating plant. Given the large local sources of emissions, the heavily populated areas to the south and west, and the potential for local stagnation of airflow, there is obviously great concern regarding air quality in the Port Moody Basin. Ozone concentrations at station T l l (Abbotsford Airport) have also been found to exceed air quality standards (CSC, 1985). In contrast to station T9, station T11 is located near a lightly-used airport in an area of largely agricultural land-use where no large local sources of emissions exist. Station T 11 is directly downwind of the major emissions source area (i.e. the greater Vancouver region) when O3 concentrations throughout the valley are highest (CSC, 1985). Hence, it is believed that the source of the photochemical pollutants that reach station T11 is the greater Vancouver region. 3. D E S C R I P T I O N
305
In order to reproduce the seasonal variation of Oa concentrations, Horowitz and Barakat fitted a second-order polynomial (regressed on Julian day number t and on t 2) to data that had first been logtransformed. The use of a polynomial to describe annual variability appears appropriate in as much as the ordinary least squares procedure objectively determines the timing of the summer maximum. Horowitz and Barakat (1979) applied a simple log transformation to the data before determining the least squares coefficients of the polynomial; here, the data are transformed via the three-parameter lognormal procedure outline by Ott and Mage (1976). A second-order polynomial is sufficient to remove any trend in the data from stations T9 and T l l . After removal of the polynomial, a highly autocorrelated sequence of data (e,) remains: et =
l n ( x , - K ) - (b o + bit
Three forecast models are described and developed using data from stations T9 and T l l during the months of May to September, 1978-1985. Only high 03 months are chosen to avoid mixing data that have different characteristics (Chock, 1986). The comparative ability of the models to forecast daily maximum 1-h average 03 concentrations is presented in section 4.
et = dplet- 1 + at
(1)
(2)
where ~b1 is a first order autoregressive coefficient, e,_ is the value of the residual series on the previous day, and a t is a normally and independently distributed sequence with a mean value of approximately zero. Hence, all trend and serial correlation have been removed and the forecast O a concentration, ~t, is described by
3.1. A deterministic/stochastic (D/S) model Horowitz and Barakat (1979) have identified the general form of a non-stationary stochastic process that is able to simulate observed pollutant concentrations without assuming independence or first-order stationarity. They examined one year of 03 data from the St. Louis Regional Air Pollution Study to demonstrate the usefulness of such a model.
~t=exp[(bo+blt+b2t2)+dPlet_l]+K.
(3)
Parameter estimates for stations T9 and T l l using data from 1978 to 1985 are presented in Table 1 (1986 data were reserved to evaluate the model's performance). Previously, this model has been used to
Table 1. Parameter estimates (1) Deterministic/stochastic model Station bo b1 T9 3.57 8.58 x 10-2 (0.02)* (2 x 10-4) Tll 3.81 4.45 x 10-2 (0.01) (2 x 10 -4)
K
b2
-2.45 × 10 -4 (6 x 10- ~) - 1.33 x 10 -4 (5 x 10- 7)
b2
0.309 (0.04) 0.433 (0.04)
*Numbers in parentheses are standard errors of estimates.
AE(B) 24:2-[
b2 t2)
where Zt is the ozone concentration on day number t, K is a constant, and b o, b 1 a n d b 2 are ordinary least squares regression coefficients. Inspection of autocorrelation and partial autocorrelation functions for the residual series et clearly suggests a first-order autoregressive equation (for all years):
OF MODELS
(2) ARIMA (1, 1, 0) Model Station ~bl T9 0.739 (0.08)* T11 0.756 (0.08) (3) TEMPER Model Station bo b1 T9 - 16.3 2.48 (3.7)* (0.19) Tll 1.51 1.15 (2.8) (0.14)
+
0.526 (0.02) 0.592 (O.O2)
-25 -25
306
S.M. ROBESO~qand D. G. STEV~q
forecast probabilistic aspects of Oa time series (Robeson and Steyn, 1989). The model is termed deterministic/stochastic (D/S) in reference to the seasonal polynomial (deterministic component) and the autoregressive term (stochastic component).
Mean values of estimated model parameters for stations T9 and T l l are given in Table 1. The ability of the ARIMA (1,1,0) model to forecast daily maxima of O a is discussed in section 4. 3.3. A temperature and persistence based model
3.2. An A R I M A model The time series method of Box and Jenkins (1976) is a commonly used technique for forecasting O3 concentrations (e.g. Merz et al., 1972; Sandberg et al., 1974; Chock et al., 1975; McCollister and Wilson, 1975; Simpson and Layton, 1983). A purely stochastic autoregressive integrated moving average (ARIMA) model may be developed. When given an original series L, the general form of the model is represented by w , = 471 w t - 1 + "
" " + dpp w t _ t, + a t - 0 1 a t - 1
.....
Oqar_ q
(4) where w, is the differenced original series, w, = VZ,.
(5)
Box and Jenkins (1976) define the backward difference operator as V~(,= X,- Z,- 1
(6)
which may be successively applied to the series d times. A troublesome property of the differencing filter is its 'noise amplification' (Bennett, 1974). Although in practice d rarely exceeds 2 (Box and Jenkins, 1976), differencing should be used carefully since the magnification of high frequency components may be consequential. The autoregressive (qb) and moving average (0) model parameters are determined either by maximum likelihood or least squares estimation. The effects of all factors other than past time series values are incorporated into the 'random shock' term, at. Proper model identification should result in an a, series that is normally and independently distributed with a mean value of approximately zero. Inspection of autocorrelation and partial autocorrelation functions for the 0 3 series suggests that the most appropriate model (i.e. the one with relatively small white noise variance and a small number of parameters) is an ARIMA (1,1,0) model, where the notation ARIMA (p,d,q) refers to the number of autoregressive parameters, p, moving average parameters, q, as well as the number of differencing operations performed, d. Various versions of Equation (4) were tested to find a model that accurately describes seasonal variability and maintains some degree of parsimony; however, with increasingly complex models, very small decreases in residual variance resulted. An ARIMA (1,1,0) model has the form w t + q91 w , _ x = at.
(7)
When transformed to the original series, this becomes ~t = (~b 1 + 1 ) Z , _ I - 4~1Z,- 2 + at.
(8)
The two models previously described are variations of the 'pure time series' approach, i.e. a univariate series is decomposed into various elements. In this section, a model that relates daily maximum 0 3 to 'external' variables as well as to previous O3 values is developed. Forecast models often rely heavily on regression techniques for both model formulation and estimation (Aron and Aron, 1978; Wolff and Lioy, 1978; Aron, 1980; Prior et al., 1981). Numerous independent variables commonly are entered into a stepwise regression to identify those having 'significant' partial correlations with daily maximum 0 3. Such models typically explain a non-trivial proportion of the variance in daily O a maxima; however, practical and methodological difficulties preclude their use here. First, many of the independent variables commonly used are derived from upper air data (e.g. Aron and Aron, 1978). Such data were not available for this study since no upper air stations are located in the lower Fraser Valley. Second, simply entering a large number of independent variables into a stepwise regression procedure is neither physically purposeful not intuitively appealing. When selecting independent variables to forecast 0 3 variability, factors that influence both photochemical production and atmospheric accumulation should be considered. Rudimentary examination of Oa photochemistry indicates the importance of NO, NO 2 and certain species of HCs (Fishman and Crutzen, 1978). Ozone accumulation is related to dispersive properties of the atmosphere and hence to wind speed, inversion height and stability. Solar radiation (in particular, u.v. wavelengths) also affects O3 variability in a well-known fashion. Munn (1975), however, has shown that the strength of solar radiation input across Canada is not a limiting factor in photochemical reactions during summer. Of all meteorological variables, air temperature has the strongest correlation with 03 concentrations for two reasons: (1) high air temperatures are an excellent indication of environmental conditions conducive to 03 production and accumulation (i.e. anticyclonic conditions with associated clear skies and light winds) and (2) the rate 'constants' of photochemical reactions are highly temperature dependent (OECD, 1979). As a result (at least for the lower Fraser Valley), air temperature is a reasonable surrogate for the combined effects of wind speed, wind direction, inversion height and photochemical reaction rate. Consequently, correlation and partial correlation analyses were performed to identify important independent variables from a set including NO, NO2, wind speed and its inverse, u.v. radiation,
Statistical forecast models for ozone air temperature and the previous day's O 3 concentration. Speciated HC data were not available. A combined linear relationship between daily maximum 03 concentration and (I) daily maximum air temperature (Tt) and (2) the previous day's daily maximum 03 concentration (Xt-l) explains a large proportion of the variance in Xt. This may be expressed as
307
regression of forecast variable on observed variable, mean absolute and root mean squared errors (MAE, RMSE), and systematic and unsystematic components of RMSE: [-1 n
^
271/2
RMSE~=|- ~ (P,-O,) / Lni=l
(10)
J
[-1 ~ 2]1/2 RMSEu = Lni~I(PI- Pi) ] , .
f(, = bo + bl Tt + b2z,- ~
(9)
where bo, bl and b2 are coefficients from an ordinary least squares regression (see Table 1). Explained variance (r 2) at station T9 is 0.52 while for station T l l , r 2 = 0.41. These values of explained variance are slightly lower than those obtained by Prior et al. (1981) via a similar approach; however, unlike that study, the data used here (1) contained multi-year variability, (2) did not include early morning values, and (3) were not spatially averaged. As a result, forecasts based on Equation (9) are both timely and meaningful for operational purposes. Daily maximum air temperature data were taken from 1-h average air temperature values at station T9 from 1982 to 1985 since temperature data were not available at station T9 prior to 1982. If used operationally in a forecasting mode, Equation (9) requires accurate forecasts of daily maximum temperature. The only temperature forecast verification performed in the lower Fraser Valley occurs at Vancouver International Airport where forecasts for the months of June-August of 1986 had a mean absolute error of 1.27°C (Hammond, pers. comm., 1987). Further inland (i.e. away from the moderating influence of the Strait of Georgia), temperature deviations are likely to be larger; nonetheless, daily maximum air temperature forecast errors of I°C result in 03 forecast errors of only 2.5 ppb at station T9 and 1.1 ppb at station Tll--values small enough to warrant using forecast air temperature. This bivariate regression model, Equation (9), is termed TEMPER due to its reliance of Temperature and Persistence. Its ability to forecast daily maxima of 03 is examined in the following section. 4. M O D E L PERFORMANCE
Three statistical models of the temporal variability of 0 3 were described in the previous section; however, no attempt was made to compare the relative ability of each model to forecast daily maximum O a concentrations. In this section, the forecast performance of each model is evaluated via graphical and statistical comparison with data from the year 1986--data not used in model development. In addition, all models are compared with a persistence forecast (i.e. f~t= ;(t- ~), a benchmark comparison in nearly all of the atmospheric and environmental sciences. As recommended by Willmott (1981) and Willmott et al. (1985), an array of model evaluation statistics are presented: observed and forecast means and standard deviations, intercept (a) and slope (b) of a least squares
(11)
where Oi and Pi are the observed and predicted values, respectively and Pi = a+ boy The simple linear relationship RMSE 2 = R M S E 2 + R M S E 2 indicates the very useful decomposition of the total error into systematic and unsystematic elements. Also reported is the Index of Agreement (Willmott, 1981), d, defined as d= 1
X~=~(Pi- O~)~ , , , 2; ~ : 1(lPil + IO,I)
(12)
where P [ = P i - O and O ' i = O i - O , with 0 being the observed mean value. Being dimensionless and having the limits of 0.0 (indicating no agreement) and 1.0 (indicating perfect agreement), d may be viewed as a standardized (by the variability in the predictions and observations about the observed mean) measure of the mean squared error. To forecast daily maximum O a concentrations with the D/S and ARIMA models, the random term, at, is simply set equal to its expected value: zero. Forecasts using the TEMPER model employed actual, rather than forecast daily maximum air temperature values. To test the sensitivity of TEMPER to air temperature forecast errors, a series of normally distributed random numbers (with mean of 0.0 and standard deviation of 2.0) is added to observed temperatures. Air temperatures containing these random errors are then used in TEMPER to forecast daily maximum 03 concentrations. Scatterplots of observed vs forecast daily maximum 03 concentrations are presented in Fig. 2 (a)-(d) for station T9 and Fig. 3 (a)-(d) for station T11. Scatterplots for the TEMPER with error forecasts are not included since they are nearly indistinguishable from those of the TEMPER forecasts: Model evaluation statistics are summarized in Table 2 for comparative purposes. An examination of Table 2 reveals a few noteworthy facts. First, the persistence forecast was better than that of the D/S model and nearly the same as that of the ARIMA model at both stations T9 and Tll. Second, for all models, the slope and intercept of the regression line indicate that low concentrations are being overpredicted and high concentrations are being underpredicted; apparently, none of the models are sufficiently dynamic to capture the rapid variability in the 03 time series. Although nearly all statistical models underestimate variance (since they are usually some form of non-magnifying filter through data), the D/S model performs particularly poorly in this regard. Third, the random noise added to the daily maximum
308
S.M. ROBESONand D. G. STEYN (o) o/s
(b)
140.0
140.0.
O. 120.0 O. v
Q. 120.0,
o
O.
C 100.00
(~ 100.0.
80.0-
o
C
o
.o•
oo
,
C 0
60.0-
o*" o o o
.
o
o
o
,o0 o
,
o
o
O .~L C 0 ~)
o
o .o
o
o
20.0' o
o.0 0.0 " 20:0
40:0 "60:0
"80:0
0.0
" I0().0" 12().0" 140.0
20.0
0.0
Observed
Concentrotion
(ppb)
(d)
140.0-
o
80.0-
o
o
,"
d)oO o Y
"°°°
o o 6, %0 . ,
60.0-
0
C)
120.0
140,0
(ppb)
o
o
oo
o
100.0.
o o o
O C 0 U )
o
°°° o
o
•o.o
o oo
o
ooo
i
20.0
-
,
40.0
Observed
-
,
60.0
•
1
-
80.0
Concentrotion
,
100.0
-
,
120.0
-
140.0
(ppb)
o
oo
°o°
20.0. o
oo
~°° °o i
o
o o o
~
o
O U..
o.6
60.0.
O
o
. o
80.0
0 U
0,0
100.0
CL 120.0 CL
o /
c- 100.00
0.0
80.0
Persistence
z"
Q. 120.0O.
20.0-
60.0
Concentrotion
140.0 ,
z"
~ 0 b-
40.0
Observed
(c) T E M P E R
40.0 -
-'" ..o" ° °
40.0"
O LL
"~ 0
o
o o o ° o oo o
O U
0
r-
. e
o o
;o
60.0-
,
o
o
60.0.
O
Z:oO
"~)'oo
20.o O o ~ .
.-.-
b_
C
ARIMA
°° o °°°
o
0.0 0.0
" 2010
" 4010
Observed
"60:0
"8010
Concentration
1 0 0 . 0 " 12().0
140.0
(ppb)
Fig. 2. Scatterplots of forecast vs observed daily maximum ozone concentrations for the summer of 1986 at station T9 using: (a) D/S, (b) ARIMA (1,1,0), (c) TEMPER and (d) persistence models. Solid line indicates perfect prediction while the dashed line is an ordinary least squares regression estimate.
temperature data has had little effect on the forecast performance of the TEMPER model. A traditional time series approach (decomposition of the variability into trend, cycle and stochastic elements) was taken to develop the D/S model. However, the D/S model did not perform well--the naive persistence forecast is clearly better with the data used here. The D/S model contains the most linear systematic error as indicated by its RMSEs; the ratio MSEJMSE ~ 0.70 indicates that 70% of the total error is systematic. It is likely that appropriate modifications (e.g. relating the 'whitened' series, a,, to a 'whitened' air temperature series) would improve its forecasting ability. However, by explicitly including a polynomial representation of the seasonal cycle of 03,
high and low values in the 03 series are modulated. This is accomplished when the polynomial is subtracted from the current 0 3 concentration and the residual (rather than the entire value as in the ARIMA, TEMPER, and persistence forecasts) is used as the dynamic component in the forecast. These results also indicate that the utility of the D/S model as a conditional probability density function (Robeson and Steyn, 1989) may be compromized. In other disciplines such as economics where time series analysis also plays an important role, the traditional approach has been found to be inadequate in many situations (Granger, 1980). Hence, when using univariate methods, alternative strategies such as adaptive parameter estimation may be necessary to empirically model
Statistical forecast models for ozone (o) o/s 100.0 O. (3.
80.0
C 0
cO
60.0-
0
(..)
o o
~o
.-"
omo o
o
C 0
o
o
Qo
o o
U tO
o
0
°. . ,
0 b_ 0.0 0.0
0 (.)
2o.o
40.0
6o0
ao.o
0.0 0.0
~oo.o
0
20.0
(d)
TEMPER
o
"'"
oO ooo
~.o'o
40.0
eo.o
Bo.o
100.o
Persistence
100.0
I
°° ~
80.0.
80.0-
°
oo
o o 60.0-
.~
(9
60.0
g
o
0
40.0
4O.0.
20.0
20.0
.'"
o
. o.,
o
o
o o
o °
o ~ o ~ o . "1~,~i~ o°
~
o o
o o oo o , °° ~ , ~ o ° o
"~' J ' ~ ' C
O
o
oO
0
O O
oo
Observed Concentration (ppb)
C E
o° o o
.-
2O,0
C 0 "~
, ." o . -
0 L
100.0
O. CL
o
Oo o o o L ~ , / ~c oo oo o - oOO -" oo~, ~ o -
ObseFved ConcentFation {ppb)
(c)
o
o
40,0
o
20.0-
oo
60,0
o
C)
o . -~
o o
40.0- -o o
o
80,0
o
o
o
O.
o
o .-'"
ARIMA
100.0-
-'°"
o o ®o
o o
oo o
o
C
(J c
o
(b)
309
To
o -o
o
o
0
0.0 0.0
20:0 Observed
' 40.0
" 60.0 '
Concentrotion
" 80.0 '
100.0
(ppb)
0.0 0.0
2olo 4 o l o Observed
~o:o
Concentration
6olo
~oo.o
(ppb)
Fig. 3. Scatterplots of forecast vs observed daily maximum ozone concentrations for the summer of 1986 at station T11 using: (a) D/S, (b) ARIMA (1, 1,0), (c) TEMPER and (d) persistence models. Solid line indicates perfect prediction while the dashed line is an ordinary least squares regression estimate.
systems (such as the photochemical one) that vary on short timescales. To gain an understanding of why the ARIMA and persistence models performed similarly, it is useful to view the two models as filters: the original time series is transformed into an uncorrelated white noise sequence. The persistence model transforms a time series in exactly the same way as the first-order differencing filter described in section 3.2 does, while the ARIMA (1,1,0) model is simply an optimally chosen (in this case by maximum likelihood estimation) weighted average of previous values in the Oa series. Comparison of the amplitude spectra of the ARIMA(1, 1,0) and persistence filters (Fig. 4) reveals the reason for the similar forecasting ability--their amplitudes in the frequency domain are nearly identical (as one might
also suspect based on the magnitude of the autoregressive coefficients). Phase spectra of the two filters are similar as well. The widespread use of the Box-Jenkins (ARIMA) method in the atmospheric and environmental sciences (McLeod et al., 1977; Bennett, 1979) results from its ability to forecast state variables while maintaining model simplicity (the models presented here have only one parameter, ~bt, to estimate). However, the performance of such models should always be compared--from the perspective of model evaluation, not simply in terms of the residual variance in the fit--with that of persistence. Regardless of simulated forecast errors in daily maximum air temperatures, the TEMPER model appears to be superior to the other models as indicated by the array of model evaluation statistics (in
S. M. ROBESON and D. G. STEYN
310
Table 2. Model evaluation statistics Station T9
(1)
Observed mean Observed standard deviation Predicted mean Predicted standard deviation Regression line intercept slope Mean absolute error Root mean square error systematic component unsystematic component Index of agreement
Station T l l Observed mean Observed standard deviation Predicted mean Predicted standard deviation Regression line intercept slope Mean absolute error Root mean square error systematic component unsystemetic component Index of agreement
Model* (3)
(2) - -
41.3
(4)
(5)
48.1 18.3
41.9 21.0
21.9 0.632 12.9 16.2 10.2 12.6 0.821
18.7 0.560 13.2 19.7 9.3 17.4 0.743
- -
--21.1-49.0 17.5
42.7 11.4
41.7 19.3
29.7 0.313 13.2 17.3 14.5 9.3 0.678
21.0 0.500 13.6 19.3 10.5 16.2 0.733
(1)
(2)
Model (3)
(4)
(5)
39.6 9.4
41.5 15.2
41.3---16.6~ 43.6 11.3
43.0 11.6
41.5 16.5
26.0 0.328 9.8 13.6 11.3 7.6 0.686
21.5 0.484 10.7 15.5 8.5 12.9 0.722
23.3 0.476 9.7 12.2 8.8 8.5 0.792
17.9 0.570 10.4 15.3 7,1 13.5 0.752
22.7 0.636 12.3 15.6 10.8 11.2 0.832
23.8 0.480 9.3 12.0 8.9 8.1 0.798
*(1) D/S, (2) ARIMA, (3) TEMPER, (4) TEMPER with error, (5) persistence.
2.0
PURE PERSISTENCE
1.5 t-,.J
1.0
0.5
0.0 0'
;o
40 100 120 140 ' 2o 8'0 ' ' ' LINEAR FREQUENCY (CYCLES PER YEAR)
160 '
180 '
Fig. 4. Amplitude spectra of the filtering process represented by pure persistence and the ARIMA (1,1,0) models at stations T9 and Tll. particular, note the low M A E a n d high Index of Agreement). P e r f o r m a n c e (and the original fit) is s o m e w h a t p o o r e r at station T11 a l t h o u g h this simply m a y be a result of the t e m p e r a t u r e d a t a being d r a w n
from station T9. T E M P E R also forecasts reasonably well on days with high 0 3 concentrations, a n important prerequisite for a n operational forecast model. Since T E M P E R is a regression model, it requires
Statistical forecast models for ozone (1) accurate forecasts of predictor variables (i.e. daily maximum air temperature) and (2) periodic reestimation of model parameters as new data become available. All of the models appear to have some difficulty forecasting high 0 3 values; however, it should be noted that none of these models were developed to explicitly handle extreme values. Extreme value prediction is often handled in some ad hoc fashion (Prior et al., 198 I)---no such corrections were attempted here. Given the hazards of exceeding epidemiologically determined standards, though, the models' performance with high concentrations suggests that a different form of model (e.g. quadratic or multiplicative rather than a linear dependence) is worth exploring. 5. CONCLUSIONS Three statistical models and persistence have been used to forecast daily maximum 0 3 concentrations in the lower Fraser Valley of British Columbia. Through graphical and statistical comparison with data not used in model development, the relative utility of each model has been demonstrated. Of the three statistical models, a regression model (TEMPER) based on daily maximum air temperature and the previous day's 0 3 concentration provides the most accurate forecast. Incorporating physically meaningful variables such as air temperature clearly gives a regression model an advantage over univariate methods. An A R I M A model forecasts no better than pure persistence while a traditional time series decomposition (the D/S model) performs the worst with these data. Hence, given reasonably accurate forecasts for daily maximum air temperature, it is suggested that T E M PER be used for forecasts of daily maximum 0 3 concentrations in the lower Fraser Valley. The best 'rule of t h u m b ' - - a n d certainly the simplest--forecast is persistence. The models do appear to have some difficulty forecasting very high concentrations. Although ad hoc corrections may improve local forecasts, such models will likely be site specific. It is suggested that different forms of models (e.g. including quadratic or multiplicative dependencies) be developed and evaluated for general applicability. Acknowled#ements--The Pollution Control Section of the Greater Vancouver Regional District kindly provided the data used in this study. Thoughtful comments from Cort Willmott (University of Delaware) and two anonymous reviewers are appreciated. This research was supported by a grant from the Natural Science and Engineering Research Council of Canada and by a University of British Columbia Graduate Fellowship to SMR. REFERENCES
Aron R. (1980) Forecasting high level oxidant concentrations in the Los Angeles basin. J. Air Pollut. Control Ass. 30, 1227-1228. Aron R. and Aron I.-M. (1978) Statistical forecasting models:
311
II. Oxidant concentrations in the Los Angeles basin. J. Air Pollut. Control Ass. 28, 684-688. Bates D. V., Bell G. M., Burnham C. D., Hazucha M, Mantha J., Pengeily L. D. and Silverman F. (1972) Shortterm effects of ozone on the lung. J. appl. Physiol. 32, 176-181. Bennett R. J. (1974) Process identification for time series modelling in urban and regional planning. Regional Stud. 8, 157-174. Bennett R. J. (1979) Spatial Time Series. Pion, London. Box G. E. P. and Jenkins G. M. (1976) Time Series Analysis, Forecasting and Control. Holden-Day, San Francisco. Chock D. P. (1986) Statistics of extreme values of air quality--a simulation study. Atmospheric Environment 10, 1713-1724. Chock D. P., Terrell T. R. and Levitt S. B. (1975) Time-series analysis of Riverside, California air quality data. Atmospheric Environment 9, 978-989. Concord Scientific Corporation (CSC) (1982) Vancouver Oxidant Study, Air Quality Analysis 1978-1981, Prepared for Environment Canada, Environmental Protection Service. Concord Scientific Corporation (CSC) (1985) Vancouver Oxidant Study, Air Quality Analysis Update 1982-1984. Prepared for Environment Canada, Environmental Protection Service. Dickson D. R. and Quickert N. (1975) The chemical composition of photochemical air pollution. In Photochemical Air Pollution: Formation, Transport and Effects. National Research Council of Council, Publication No. NRCC 14096 of the Environmental Secretariat. Fishman J. and Crutzen P. J. (1978) The origin of ozone in the troposphere. Nature 274, 855-858. Goldsmith J. R. and Nadel J. A. (1969) Experimental exposure of human subjects to ozone. J. Air Pollut. Control Ass. 19, 329-330. Granger C. W. J. (1980) Forecasting in Business and Economics. Academic Press, New York. Heck W. W., Adams R. M., Cure W. W., Heagle A. S., Heggestad H. E., Kohut R. J., Kress L. W., Rawlings J. O. and Taylor O. C. (1983) A reassessment of crop loss from ozone. Envir. Sci. Technol. 572A-581A. Horowitz J. and Barakat S. (1979) Statistical analysis of the maximum concentration of an air pollutant: effects of autocorrelation and non-stationarity. Atmospheric Environment 13, 811-818. McCollister G. M. and Wilson K. R. (1975) Linear stochastic models for forecasting daily maxima and hourly concentrations of air pollutants. Atmospheric Environment 9, 417-423. McLeod A. I., Hipel K. W. and Lennox W. C. (1977) Advances in Box-Jenkins modelling 2. Applications. Water Resour. Res. 13, 577-586. Merz P. H., Painter L. J. and Ryason P. R. (1972) Aerometric data analysis--time series analysis and forecast and an atmospheric smog diagram. Atmospheric Environment 6, 319-342. Munn R. E. (1975) The oxidant climatology of Canada. In Photochemical Air Pollution: Formation, Transport, and Effects. National Research Council of Canada, Publ. No. NRCC 14096 of the Environmental Secretariat. Mustafa M. G. and Lee S. D. (1976) Pulmonary biochemical alterations resulting from ozone exposure. Ann. Occup. Hyg. 19, 17-26. Organization for Economic Co-operation and Development (OECD) (1979) Photochemical Oxidants and Their Precursors in the Atmosphere. OECD, Paris. Ott W. R. and Mage D. T. (1976) A general purpose univariate probability model for environmental data analysis. Comp. Oper. Res. 3, 209-216. Prior E. J., Schiess J. R. and McDougal D. S. (1981) Approach to forecasting daily maximum ozone levels in St. Louis. Envir. Sci. Technol. 15, 430-436.
312
S.M. ROBESON and D. G. STEYN
Reich P. B. and Amundson R. G. (1985) Ambient levels of ozone reduce net photosynthesis in tree and crop species. Science 230, 566-570. Robeson S. M. (1987) Time series analysis of surface layer ozone in the lower Fraser Valley of British Columbia. M.Sc. Thesis, The University of British Columbia, Vancouver, British Columbia, Canada. Robeson S. M. and Steyn D. G. (1989) A conditional probability density function for forecasting ozone air quality data. Atmospheric Environment 23, 689-692. Sandberg J. S., Robinson L. H. and DeMandel R. E. (1974) Box-Jenkins stochastic modeling in regional forecasting of photochemical oxidant. In AMS Conference on Atmospheric Turbulence and Diffusion, 1974. American Meteorological Society, Boston.
Simpson R. W. and Layton A. P. (1983) Forecasting peak ozone levels. Atmospheric Environment 17, 1649-1654. Steyn D. G. and Oke T. R. (1982) The depth of the daytime mixed layer two coastal sites: a model and its validation. Boundary-Layer Met. 24, 161-180. Willmott C. J. (1981) On the validation of models. Phys. Geog. 2, 184-194. Willmott C. J., Ackleson S. G., Davis R. E., Feddema J. J., Klink K. M., Legates D. R., O'Donnell J. and Rowe C. M. (1985) Statistics for the evaluation and comparison of models. J. geophys. Res., 90, 8995-9005. Wolff G. T. and Lioy P. J. (1978) An empirical model for forecasting maximum daily ozone levels in the northeastern U.S.J. Air Pollut. Control Ass. 28, 1034-1038.
