Abstract ll; 1. Introduction;; ll; 2 Methodology;; ll; 3 Data;; l (a) Rainfall;; l; (b) Sea Surface Temperature ;; ll; 4 Principal Component Analysis ;; l; (a) Rainfall;; l (b) Sea Surface Temperature; ll; 5 Preliminary Analysis;; ll; 6 Forecast Model;; l; (a) Single predictor model;; l; (b) Two potential predictors;; l; (c) Full SST data set - multiple predictors.;; l; (d) Independent Forecasts;; ll; 7 Conclusions. ;
References;; ll; Appendix A;; ll; Appendix B;;
Extensive testing of sea surface temperature anomaly (SSTa) patterns as predictors of Australian seasonal rainfall variations is carried out. Ten rotated principal components covering the Indian and Pacific Oceans are used as predictors. Individual monthly values are used with a one and three month lead time, eg November and January SSTa are used to forecast March to May rainfall. Current values of the principal component amplitudes are calculated by projecting the BMRC/NMC or NCEP Washington SST analysis onto the set of SST principal components. ;
Hindcast experiments with two and five degree area averaged rainfall suggested that skill is improved by forecasting on larger scales. Nine large-scale regions covering the Australian continent are defined by a rotated principal component analysis of a grided one degree rainfall dataset. Hindcasts are then made of the amplitudes of these nine principal components, and the principal component loadings are used as weights to project these back to the original one degree grid points. ;
Assessment of hindcast skill has been carried out using the cross validation technique. This involves a series of forecasts over all the available data, leaving one (or more) observations out each time. Hindcasts experiments, using linear discriminant analysis, are carried out on a variety of predictors, ranging from the Southern Oscillation Index (SOI) alone, the SOI and an Indian Ocean Index together, the first two SST principal components and finally all ten SST components at two lagged times. In all cases the models are restricted to the best combination of two predictors, with selection of this best combination of predictors done using a cross model validation procedure. This methodology produces much lower, but more realistic, estimates of skill than one based on the dependent data. Linear Error in Probability Space (LEPS) skill scores is used to extensively in both the predictor selection and model validation phases. ;
Experimental real-time forecasts over the period from January-March 1994 to June-August 1997 indicate improved skill over parts of southern Australia, compared to forecasts using the Southern Oscillation Index alone. ;
Past research has highlighted the significant impact of climate variations on the Australian economy, and the potential value of timely forecasts of these climate variations. Current predictions, based on the state of the El Nino / Southern Oscillation (ENSO) show some skill (above chance, persistence or climatology) of forecasts of seasonal climate variations.
To improve these forecasts a number of projects are being conducted in the Bureau of Meteorology Research Centre. This report documents current empirical work aimed at utilising global scale Sea Surface Temperature anomaly (SSTa) patterns as predictors of inter-annual Australian climate variability, particularly seasonal rainfall. Related research aims to document and understand the intraseasonal variability of Australian rainfall and the impact of this variability on seasonal prediction, and test the ability of General Circulation Models forced by predicted SST anomalies to forecast seasonal rainfall variations.
The prediction of rainfall (or some other climatological parameter) by statistical or empirical methods is feasible if there is a lagged relationship between the rainfall and a suitable predictor. Such a relationship exists in some seasons between eastern Australian rainfall and the Southern Oscillation (SO; McBride and Nicholls, 1983). This relationship provides the basis for a Seasonal Climate Outlook (SCO) issued by the National Climate Centre (NCC) of the Australian Bureau of Meteorology. The major emphases of this service, which uses the Southern Oscillation Index (SOI; Troup, 1967) are seasonal (3 monthly) rainfall forecasts for all of Australia. Forecasts are also made of the timing of the onset of the north Australian wet season (monsoon) and the expected number of tropical cyclones, which are also related to the SOI.
Since the commencement of the SCO, various forecast techniques have been employed to extract maximum information from the SOI. Initially actual seasonal rainfall was predicted using simple linear regression with the past season (three months) value of the SOI as the predictor. This predicted rainfall was then expressed as a decile value. These categorical forecasts were replaced by probabilistic forecasts in 1992, partly in response to user's requests for some error or probability estimates to be included with the forecast. The forecasts were produced using linear discriminant analysis with one of two possible predictors; the past three month seasons value of the SOI, or the three month trend in the SOI.
Currently the outlooks are based on an optimal combination (Casey, 1995) of climatology, persistence and phase categories of the SOI (Zhang and Casey, 1992). In addition further qualitative guidance is provided by an alternative SOI phase technique (Stone, 1992) and an analogue selection technique (Drosdowsky, 1994).
The relationships between northern and eastern Australian rainfall and the SO are weak, however, during the late southern summer and autumn period, when ENSO events tend to decay. Over the western third of the continent relationships are weak throughout the year, except for the southwest corner in autumn. These limitations of the SO have lead, in recent years, to attempts to find other independent predictors. Nicholls (1989) and Drosdowsky (1993a,b) have documented a major mode of winter rainfall variability, consisting of a broad band extending across the continent from the northwest of Western Australia, through central Australia to the south coast. Nicholls (1989) showed that this pattern was related to the gradient in sea surface temperature (SST) between the central Indian Ocean and Indonesia. This anomaly pattern has been called the Indian Ocean Dipole. Drosdowsky (1993c) showed that the strongest precursor to this dipole pattern was found in the SSTa in the eastern Indian Ocean near the west coast of Australia during the preceding summer and autumn.
Large scale SSTa patterns have been shown to be related to seasonal rainfall anomalies in other regions of the globe. Ward and Folland (1991) have linked rainfall anomalies in the Sahel and northeast Brazil with large scale patterns of Pacific and Atlantic Ocean SST anomalies. Barnston (1994) has used Canonical Correlation Analysis (CCA) with near global SST patterns as well as Northern Hemisphere 700 hPa geopotential height fields as predictors of North American and European precipitation and surface temperature.
Development of an empirical forecasting scheme involves a number of essential steps.
(i) identification of predictands and possible predictors
(ii) construction of the statistical model, including predictor selection when necessary, and
(iii) verification or estimation of model skill.
Various techniques have been used to find coupled patterns in climate data (Bretherton et al, 1992; Barnston, 1994). Previous studies of the relationships between Australian rainfall and near-global SSTa (Nicholls, 1989; Drosdowsky, 1993a,b) had reduced the dimensionality of the rainfall data using principal component analysis, and searched for teleconnections using simultaneous and lagged linear correlations between the principal component time series and the SST. From these teleconnections a number of simple indices of SSTa averaged over large areas similar to the "El Nino" indices (Reynolds, 1982) have been identified as potential predictors, eg Nicholls (1984) North Australia / Indonesia Index and the Indian Ocean Index (Drosdowsky 1993c).
A more systematic approach is used here to find potential predictors in the SST data. The modes of variability in the SSTa are summarised using principal components which are then related to the rainfall. This has a number of advantages, over the combined analysis techniques, such as canonical correlation or singular value decomposition (SVD), most importantly that other predictors, eg atmospheric circulation patterns, and predictands, eg maximum or minimum temperatures can be easily introduced.
Once a model is selected its parameters are estimated from all or part of the available data. The model is then tested on an independent set of data to obtain an estimate of its skill. If all the available data is used to fit the model parameters, a valid estimate of the skill cannot be obtained. The available data must therefore be split into a dependent set used for model development and an independent or validation set. Two methods of developing a forecast model have been used in this work. Initially a traditional approach of model selection or specification was used. In this approach the data set of length N is split into two parts; a training set of (about) 2N/3 is used to develop the model and estimate any parameters, and the remaining third of the data is used as a verification set on which independent hindcasts are performed and assessed. In general, the forecast errors will be larger for this verification data set than for the development set.
If the number of data available (N) is not very large, as is usually the case in climate prediction, the estimates of skill obtained on the verification set may not be an accurate assessment of the true skill. A more accurate assessment can be obtained through the technique of cross-validation. There are two uses of the cross validation methodology (Stone, 1974; Michaelsen, 1987; Elsner and Schmertmann, 1994), that of model assessment or single cross validation as outlined in Appendix A, and that of model choice or double cross validation as described in Appendix B.
Use of the single cross validation technique for estimating the skill of a forecast model involves a series of independent forecasts over all the available data. This is achieved by leaving one (or more) observations out each time. The model parameters are then recalculated and a hindcast produced for the missing observation. If the observations are strongly auto-correlated blocks of two or more observations may need to be omitted to maintain the independence of the forecast and development data subset.
The methodology validates a particular type of model rather than the specific set of parameters. Since there are now N hindcasts, this method should give more accurate estimates of the true hindcast skill than the traditional methodology of splitting the data set into single training and verification sets and validating both the model type and the particular set of parameters.
In the double cross validation, different models may be assessed, eg selecting a subset of predictors from a larger pool of potential predictors. Firstly the type of model to be considered is defined, and the inner cross validation loop is used to find the best model of this type for a data set with one year left out. In the outer loop this best model is used to produce a hindcast for the left out case. The skill of the class of models is then evaluated over the full set of hindcasts.
Within these development and verification frameworks, two statistical methods have been examined to forecast rainfall using the SSTa; multiple linear regression (MLR) and multiple discriminant analysis (MDA). Regression forecasts of grid point rainfall amounts were verified using absolute (root mean square; RMS) errors. Probabilistic forecasts for rainfall categories (quintiles or terciles) were verified in two ways; a simple hits score in the preliminary testing described below, and the Linear Error in Probability Space (LEPS) score described by Potts et al, 1996.
(a) Rainfall
Earlier works noted above (Nicholls, 1989; Drosdowsky 1993a,c) used the Australian district rainfall data set. This has a number of known deficiencies, specifically the variable size of the districts and the changing stations within each district over time (Jones and Beard, 1998). To overcome these problems, two new rainfall data sets have been used. The first, referred to as the box averaged rainfall, was created at one degree spatial resolution and monthly temporal resolution by averaging all available data from 1910 to 1993 into one degree squares (Drosdowsky, 1993d). More recently, a new rainfall data set at one quarter degree spatial and monthly temporal resolution has been created by the National Climate Centre (Jones and Weymouth, 1997), by analysing all the available station data using a successive correction scheme. A reduced resolution subset (at one degree spatial resolution, and referred to as the grided rainfall data) has been used in this study.
The climatology and interannual variability of the box averaged and grided rainfall data sets are virtually identical, and the results of some of the initial analyses performed on the box averaged data have been assumed to be valid also for the grided rainfall data. The spatial and temporal variability, as defined by an Principal Components Analysis will be described in the following section.
(b) Sea Surface Temperature
For the initial analysis the UK Meteorological Office Global Ice and Sea Surface Temperature data set (GISST, version 1.1; Parker et al, 1995), from January 1949 to December 1991, has been used. The variablilty of this data set is examined in the following section using Principal Component Analysis. A new version of the GISST data set (version 2.1) extending to 1994 became available during the course of this work. The SST principal components were recalculated using this data but have not been utilised due to apparent problems in no-data areas due to the interpolation and reconstruction method.
To extend the SST data from December 1991 to the present, two operational analyses are used. The first of these, the National Centre for Environmental Prediction (NCEP) optimum interpolation analysis is available from 1982 to the present (Reynolds and Smith, 1994), while a locally produced analysis (BMRC/NMC Melbourne) is available from June 1993 (Smith, 1995). All these data sets have 1o by 1o spatial resolution. The 1o GISST data for the 43 years from 1949 to 1991 were reduced to a 5o by 5o degree grid, with 848 boxes covering the region from 60oN to 55oS and 30oE to 70oW. The operational analyses are also reduced to the 5o by 5o degree grid and projected onto the GISST1.1 principal components using the scoring coefficients matrix.
For both the rainfall and SST data sets standardised monthly anomalies were calculated at each gridpoint by removing the full data set monthly mean and dividing by the monthly standard deviation, as described in Drosdowsky (1993a). The resulting time series have zero mean and approximately unit variance and, hence, the correlation and covariance matrices are virtually identical. The rainfall and SST fields have been analysed using Varimax rotated S-mode principal component analysis as described in Drosdowsky (1993a). The principal component analysis was performed using IMSL subroutines.
The relationships between Australian rainfall and some simple predictors, eg the SOI (McBride and Nicholls, 1983; Drosdowsky and Williams, 1991) or Indian Ocean SST (Nicholls, 1989; Drosdowsky, 1993c), are known to vary with season and location. A forecast scheme based on SST would be expected to exhibit similar seasonal variation, so that a stratification of the data sets is necessary. The stratification of the monthly data can be done before or after the Principal Component Analysis. Drosdowsky (1993a) examined the effect of seasonal stratification on a non-overlapping seasonally averaged district rainfall data set and concluded that there was little difference in the two analyses over most of the continent. Similar tests were performed using both monthly and two or three month, non-overlapping, seasonally averaged data, for both of the rainfall data sets and the SST data. In all cases the spatial patterns (principal component loadings) produced by the monthly and seasonal data sets were virtually identical, as were the resulting principal component time series when the monthly series were averaged into the appropriate seasons.
In addition, although forecasts are currently issued for a three month season, the feasibility of one and two month forecasts should also be examined. Therefore, only the unstratified monthly principal components are used in the following analyses, since the seasonal values can be easily reproduced by averaging the appropriate monthly time series.
(a) Rainfall Principal Components
The analysis of seasonal district rainfall anomalies described in Drosdowsky (1993a) has been repeated using the box average and grided data sets, and for a slightly longer period, 1950 to 1993 (1995 for the grided data). In both cases the eigenvalue spectrum (Figure 1) suggests breaks at four and at nine components. Since we wish to retain as much of the variance as possible, the higher breakpoint was selected. When the first nine components are rotated using the Varimax procedure (Figure 2), the eight regions found earlier (S1 to S8 in Drosdowsky, 1993a) are broadly reproduced, although the relative importance is altered. Although the box average and grided data sets also have differing importance for each component, the spatial patterns remain more closely related to each other than to the district rainfall data. This is largely due to the uniform coverage of the one degree grid of these data sets compared to the rainfall districts. The additional component (Rain8 in Figure 2) in the box average and grided analyses covers an area over the interior of eastern Australia, reflecting the change in data density; there is relatively more variance in the inland area since there are more one degree resolution grid squares than rainfall districts in this region.
The strongest resemblance between all three analyses is found in the south and west of the continent. The major component over Queensland (S3) is restricted to the northern part of that state in the box average and grided data (Rain7 in Figure 2), while the minor east coast component (S8) is now the dominant component over southern Queensland and northern New South Wales (Rain1 in Figure 2). The major component in the district analysis, the southeast component (S1) has much weaker loadings and covers a much smaller area in both the box average and grided analyses (Rain9 in Figure 2).
(b) SSTa Principal Components
The SSTa analysis was performed for three different areas; the whole Indian /Pacific Ocean domain, and the two oceans separately. The spatial patterns (principal component loadings) produced by the monthly and seasonal data sets were virtually identical for all three areas and, as was the case for the rainfall data, so were the resulting principal component time series when the monthly series were averaged into the appropriate seasons.
Apart from the first component, the combined basins analysis, however, is not simply the sum of the two separate analyses. The time series of the first component in all three ocean areas is significantly correlated with the SOI, and the spatial patterns represent the mature phase of an El Nino / Southern Oscillation event.
In the combined analysis, the much greater number of grid points of the Pacific Ocean tends to dominate over the contribution from the Indian Ocean. The percentage of explained variance is very similar for the Pacific Ocean and combined basins, which is to be expected since the Pacific Ocean basin covers the majority of the combined oceans area. In each case the first principal component accounts for about 15% and the first 10 components about 42% of the total variance. For the Indian Ocean the first component explains only 11% of the variance, however, the first 10 account for 55% of the basin variance.
The principal component analysis was undertaken primarily for data reduction purposes, and the unrotated principal components have the desirable property, for predictive purposes, of being uncorrelated in time. However, this is offset by the orthogonality constraint on the spatial patterns, which makes it difficult to relate most of the higher SSTa components with complex spatial patterns to the rainfall in a physical manner. Hence rotation of the SSTa principal components was also performed to localise or regionalise the variance of each of the components. This process should make it easier to relate the SSTa and rainfall anomalies in a physical rather than purely statistical manner.
The scree plot of eigenvalue against eigenvector number (Figure 3) shows possible breaks at 5 and 12. Components 11 and 12 overlap substantially, using the errors estimates given by North et al (1982), while there is little overlap between components 12 and 13. As with the rainfall data, the larger number was selected and 12 components were retained to ensure that the variance associated with the regional and seasonal patterns, such as the Indian Ocean dipole was retained.
Principal components were rotated using the VARIMAX criterion. Spatial patterns (the principal component loadings) and time series (the principal component scores) are shown in Figure 4 for the GISST data, and referred to as SST1 through SST12. SST1 represents the mature phase of an El Nino / Southern Oscillation (ENSO) event, with largest (positive) anomalies in the central and eastern equatorial Pacific and in the Indian Ocean, and weaker negative loadings in the Pacific Ocean midlatitudes around New Zealand.
The Indonesia / Indian Ocean dipole pattern does not appear in the all-months, combined oceans analysis. This may be due to it being primarily a winter-time phenomenon, and also the Varimax rotation which tends to produce components with single localised areas of high loadings. The non-tropical areas of the Indian Ocean are covered by SST2 and SST6, the first located just west of the Australian continent and extending northwest to the central equatorial region. The time series of this component is strongly correlated with the Indian Ocean Index devised by Drosdowsky (1993c), and shown to be a precursor to the dipole pattern and strongly related to Australian early winter rainfall.
SST6 is located further west and may be affected by the Algulhas current, although the strongest temperature gradients associated with this current are found further to the southwest, at the edge of the analysis domain. Both these components (especially SST6) show a strong warming trend over the 43 year record. This is also evident in SST12 which covers the South China Sea and Indonesian sea (note that the loadings on this component are predominately negative, hence the time series shows a downward trend). SST3, SST4 and SST8 are WSW to ENE orientated bands across the North Pacific. The western end of SST3 lies over the west Pacific warm pool while SST8 extends eastward from Japan in the region of the Kuroshio current. SST5 and SST7 are located in the far north Pacific, SST5 in the Gulf of Alaska and SST7 along the North American coast.
The final three components cover the South Pacific, with SST9 extending from around northern Australia to the north of New Zealand. Components SST10 and SST11 have their largest loadings in regions with little data, especially prior to the availability of satellite SST data in 1982. These two components were not used as predictors.
5. PRELIMINARY TESTING ;
Initial hindcast experiments were performed using 42 years of data from 1950 to 1991 to determine which SST components would contribute to the forecast skill, the effects of different lag times, and rainfall data resolution. These experiments used :
(i) stepwise multiple linear regression;
(ii) the traditional verification approach with a training period of 31 years from 1950 to 1980 and an 11 year test period from 1981 to 1991;
(iii) individual grid point rainfall data at one, two and five degree resolution, to determine the effect of different resolutions; and
(iv) SOI and Indian Ocean Index as well as one or two months SSTa, with a lead time of zero (eg January-February-March data is used to predict April-May-June rainfall) or one month (eg January-February-March data used to predict May-June-July rainfall).
The expectation of most users of seasonal predictions is for forecasts expressed in probabilistic terms. Although reasonable predictions were obtained using multiple linear regression, the results were not readily applicable to operational forecasting, since they are expressed in terms of categorical values. Linear discriminant analysis was tested as a forecast method to predict probabilities of rainfall categories using the set of three best predictors, as determined from the multiple linear regression. Forecasts were performed for the two or five degree, monthly or seasonal rainfall data classified into either three (terciles) or five (quintile) categories.
The major results of these experiments were:
(i) The SSTa components improved on the skill obtained from the SOI alone. The significant components were found to vary from month to month, although all contributed at some time of the year;
(ii) The best results were obtained by using a low spatial resolution (five degree grid) rainfall data and forecasting fewer categories (terciles); and
(iii) There was only a minimal decrease in skill between the zero and one month lag forecasts using the SSTa components as predictors.
The skill achieved using discriminant analysis with the best three predictors was often below that expected, given the known relationships between some of the SSTa components (eg SST1 and SST2) and Australian rainfall. This was believed to be partly due to the use of too many predictors, and also possible due to the wrong predictors; those selected by the linear multiple regression selection may not work best with the inherently nonlinear discriminant analysis procedure.
The usefulness of these forecasts is also limited by the coarse five degree resolution at which the best results were achieved. The principal component analysis described above showed that a large proportion of the total monthly rainfall variance (up to 60%) could be accounted for by the first nine components. When rotated these formed a strongly regionalised coverage of the Australian continent. These have been normalised so that the loadings represent correlations between the principal components and the original standardised monthly anomaly time series at each grid point. The squared loadings, therefore, represent the percentage of variance each component accounts for at each grid point. These squared loadings can be used to interpolate forecasts of the component amplitudes back to standardised rainfall anomalies at the original one degree resolution grid points. This procedure should improve the overall accuracy by forecasting at a very low resolution (nine components over the country) but still provide detail down to the original one degree grid resolution.
The results from the regression forecasts and original discriminant analysis procedure suggested that there was little gained by using multiple predictors. Hence the type of model chosen would consist of a maximum of two SST predictors for each rainfall component. ;
(a) Single predictor model
To establish a benchmark skill level, hindcasts were performed and assessed in the manner described in Appendix A for the 1950-1993 period using the previous seasons SOI This zero lag method is the current operational procedure. The cross-validated LEPS skill scores are shown in Figure 5. While the SOI values can be obtained immediately at the end of a month, the time required to produce and disseminate the forecast suggests that a more useful system would be to use the SOI lagged by one month The LEPS scores for this method, shown in Figure 6, show slightly lower skill overall, but a similar spatial coverage to the no-lag case. In either case the skill level is similar to the current operational system. (Beard, personal comm).
The dominant mode of SST variability over the Indian/Pacific Ocean domain, SST1, corresponds to the El Niņo and its time series is strongly correlated with the SOI. Its skill alone as a predictor of seasonal rainfall, shown in Figure 7, is comparable to that of the SOI with no lag, and superior to the SOI with one month lag. The major differences are the improved skill over northern Australia, and the slightly reduced skill over southeast Australia during spring, particularly September-October-November. A similar result was reported by McBride and Nicholls (1983), who found a strong relationship in spring between rainfall over southeast Australia and both the SOI and Darwin pressure, but not Tahiti pressure. This suggests that southeast Australia may be more strongly influenced by the local, rather than remote, ENSO related anomalies. ;
An additional independent predictor of Australian rainfall is the Indian Ocean Dipole (Nicholls, 1989). A precursor signal to this dipole pattern, called the Indian Ocean Index, was identified by Drosdowsky (1993c). This is calculated as the average SST anomaly over the area 22oS to 34oS and 90oE to 110oE. Bi-monthly values of this index of SSTa in the eastern Indian Ocean in mid to late summer were strongly correlated with early winter (May-June and June-July) rainfall over parts of southern and eastern Australia. This result is verified in Figure 8, which shows the LEPS scores obtained by using bi-monthly values of this index at one and three month lags (ie January-February and March-April anomalies are used to forecast June-July-August rainfall). A further area of significant positive LEPS scores is found over the southwest of the continent during spring and early summer, a period not examined by Nicholls (1989) or Drosdowsky (1993a,b,c).
The second mode of SST variability, SST2, has its greatest loadings in the region identified by Drosdowsky (1993c) as the precursor to the Indian Ocean dipole and used to construct the Indian Ocean Index. The seasonal and spatial variation of skill of SST2 (Figure 9) using single monthly values lagged by one and three months (ie February and April anomalies are used to forecast June-July-August rainfall) is similar to that of the Indian Ocean Index as shown in Figure 8. Restricting the predictors to a single months lags (Figure 10) results in a lower level of skill, also consistent with the results of Drosdowsky (1993c), which showed that some of the strongest relations were found at the longer lag.
(b) Two potential predictors
Inspection of Figures 6 and 8 (or 7 and 9) suggests that using both the SOI and the Indian Ocean Index, or alternatively both SST1 and SST2, as predictors should improve the overall skill in the forecasts, since either pair are essentially uncorrelated. There are two methods in which both predictors could be used. Firstly both predictors can be used together every time, ie for each rainfall principal component for every season. This simple method may degrade the forecast skill in regions where one of the predictors has little skill. This can be seen in Figure 11, which shows the skill for this method. Any degradation of skill appears to be slight, and is compensated for by increases in other areas.
An alternative method is to only use that predictor which shows skill for a particular principal component and season, including the possibility of no predictor, ie forecasting climatology. Using; Figures 6 and 8 (or 7 and 9) to (subjectively) select which predictor to use in which region leads to the skill shown in Figure 12. The main effect of this method is to reduce the areas of negative LEPS scores, and especially the magnitude of the negative scores (This is not apparent in the plots which do not contour the negative values.) This is achieved mainly through the use of climatology (ie no SST) as a predictor in regions where neither of the two available predictors shows skill, since climatology should have a LEPS skill score of zero. This improvement in hindcast skill is artificial, however, since all the available data have already been used to produce the individual predictor skill analyses, ie this skill is "in sample" or based on dependent data. If model development requires this type of selection or screening of potential predictors, then we need to perform the "double cross" validation described above and detailed in Appendix B. Note that in the case of the Indian Ocean Index in Figure 8 or SST2 in Figure 9, the simple first alternative has been used; the lagged values are treated as separate predictors and both used to forecast each rainfall principal component so that there is no predictor selection involved.
As shown in Appendix B, we begin by setting aside one year and then find the best model fit to the remaining 43 years of data. With two potential predictors, and a model allowing either or both to be used as predictors there are 4 possible combinations; climatology or no predictor, either predictor used alone or both together for each rainfall component. Since the Australia-wide rainfall pattern is represented by 9 components, there is a total of 49 different models. The search can be narrowed down considerably by noting that the rotated rainfall components are fairly strongly regionalised, hence the choice of a particular model for one rainfall component has only limited impact on the grid points most strongly represented by other components. Hence the best model fit for the entire continent is found by the exhaustive procedure of simply trying each of the 3 non-climatology predictor combinations on a particular rainfall principal component, and using climatology for all the others. Forty-three hindcasts are produced in this manner in the inner cross validation loop by leaving one year out and using the remaining 42 years in the discriminant analysis to produce a forecast for the omitted year. Having found the "best" combination of predictors for a set of data with one year removed (which may differ depending on which year has been omitted from the data set), we can now produce a prediction for that year. This procedure is repeated in the outer cross validation loop for each of the 44 years of data. The true hindcast skill for this type of model, shown in Figure 13 is then the LEPS score for these 44 hindcasts. This appears to be more similar to Figure 11, ie using both SST1 and SST2 for each prediction than to Figure 12, obtained by subjectively selecting the best predictors.
The actual selection of the best predictors using LEPS can be done in two ways. The predictands being used are the rainfall principal component amplitudes. The discriminant analysis produces probabilities for each of these which can be assessed against their observed values. However the quantity of real interest is the actual rainfall, only a fraction of which at each gridpoint is resolved by the principal components (Figure 14). When the principal component forecasts are interpolated back to the gridpoints, these can be scored against the actual rainfall. The best predictor could then be taken as that which maximises this gridpoint score in some manner, eg by averaging over the whole continent, or the number of gridpoints above some useful threshold LEPS value such as 5.0. A combination of these two criteria could also be used, with the best predictor being one that produces a positive LEPS score for the principal component, and which produces the maximum number of gridpoints with LEPS scores above 5.0. The actual method used has been to simply select the combination of predictors which give the best LEPS score for each rainfall principal component.
The final part of the model selection procedure involves fitting the best model to all the available data. This results in exactly the same set of predictors for each rainfall component and season as obtained from the subjective selection method. Since both methods use all the data to select the best predictors, the resulting skill is "in sample" and artificially inflated.
(c) Full SST data - multiple predictors
While use of two predictors (SOI and Indian Ocean SST Index, or SST1 and SST2) can improve skill over either used alone, further increases may be possible by including some or all of the other SST principal components as predictors, and allowing an extra lagged value. This increases the number of potential predictors to twenty, and the number of possible combinations to over 200, and also increases the likelihood that artificial skill will be introduced into the hindcasts, even with the double cross validation technique. The best set of up to two predictors, for each rainfall principal component and season, selected from all the development data, ie in the final stage of the double cross-validation procedure described in Appendix B, for use as predictors on new independent data are shown in Figure 15. The "in sample" or dependent data skill (Figure 16) determined by a single cross validation of this set of predictors, shows very high levels of skill for every season over almost the entire continent. This estimate is extremely optimistic, being inflated by the large number of possible predictor combinations tested. A more realistic estimate of the true hindcast LEPS skill, shown in Figure 17, is produced by applying the full double cross validation procedure outlined in the previous section and Appendix B. As with the case in the previous section, the best combination of two predictors selected from all possible predictor combinations will often differ as different years are omitted from the data set, depending on the stability of the rainfall - SST relationship.
The difference between the "in sample" dependent data (Figures 12 and 16) and the "out of sample" independent data (Figures 13 and 17) hindcast skill levels are more pronounced for the full set of twenty SST predictors than the case for just two. By increasing the total pool of potential predictors (eg, by allowing more lagged values), and the number of predictors used for each rainfall component, the "in sample" dependent data skill can be raised to (almost) any desired level. The true increase in skill, shown in the comparison of Figure 17 and Figure 13 shows the true impact on the level of hindcast skill achieved by increasing the pool of potential predictors from two to twenty. Overall there is a general increase in skill, as measured by the average LEPS score over all gridpoints, but there are also some regions where the skill is degraded when a larger number of potential predictors is used. This is caused by the large number of potential predictors, which increases the probability of selecting a predictor other than SST1 or SST2 in regions where these components have a moderate level of skill when some years are omitted. This problem can also occur in areas where there is no skill. In areas of no skill, the fitting procedure has the option of selecting climatology, but rarely does so, there almost always being at least one predictor combination from the over 200 available which gives a better fit on the dependent data in the inner cross-validation loop. In reality, there is little relationship between the predictand and the chosen predictor, and this leads to inferior "hindcasts" on the omitted case than would be the case if climatology were forecast.
This problem appears to be similar to the degeneracy noted by Barnston and van den Dool (1993) in regression forecasts. Barnston and van den Dool offer a number of possible remedies, the simplest of which would be to forecast climatology in those areas with negative skill. This procedure, however would appear to violate a basic requirement of cross-validation that the development and verification data be independent, since all the data has been examined to decide that there is no skill and climatology should be forecast.
Two alternative remedies may be to raise the selection criterion to those predictors that perform significantly better or at a "useful" level above climatology, and to omit more than one year at a time, to reduce the influence of outliers in the data. Both of these strategies have been examined, but do not appear to produce substantially different results to those shown in Figure 17.
The true skill of the SSTa as predictors of Australian rainfall lies somewhere between that shown in Figures 16 and 17, but much more likely near that shown in Figure 17. Ganeshanandam and Krzanowski (1989) performed Monte Carlo experiments on variable selection in discriminant analysis and concluded that "Method E3 (their double cross validation technique) gives consistently accurate estimation of the true error rates...", and in some circumstances ".... errs slightly on the pessimistic side".
(d) Independent Forecasts
Since no data after DJF 1993/4 have been used in the development of the forecast scheme described above, predictions for seasons after this time can be regarded as true forecasts. LEPS scores for the 42 seasonal forecasts from Jan-Feb-Mar 1994 to Jun-Jul-Aug 1997 are shown in Figure 18, aggregated over all seasons at each gridpoint (a) and aggregated over all gridpoints for each season (b). Although this is a very small sample, the seasonal cycle in skill is apparent, with highest LEPS scores during spring, particularly 1994 and 1996. Spatially, the highest skill is found over the southeast and along the northeast coast. Forecasts based on just the SOI show a similar pattern over this period, with slightly better performance on the central Queensland coast, but without the skill over the southern parts of the continent.
Overall the results show an increase in hindcast skill of the SSTa over the SOI used alone. In particular the SSTa show some skill through the so-called "predictability barrier" in late summer - early autumn. The preliminary testing showed that there is more skill in predicting on a coarser grid (smaller number of larger grid boxes), and in forecasting only three rainfall categories rather than five. This has led to the use of terciles as the predicted quantity and the use of a small number of rainfall principal components as the predictands.
Principal components are assumed to be the optimal way of representing the SSTa patterns. The best predictors from this set are then found using a cross-model validation technique (Elsner and Schmertmann, 1994), rather than the more popular combined analysis of the predictor and predictand fields together using, for example, canonical correlation or singular value decomposition. This procedure allows the use of alternatively predictors such as the simple area averaged indices such as Nino3 or the Indian Ocean Index (Drosdowsky 1993c), which have been shown to produce skilful forecasts. However the "discovery" of further predictors of this type tends to be an ad hoc procedure.
The forecasts are produced using Discriminant Analysis. This has two major desirable properties. Firstly, it can produce probability forecasts and, secondly, it is inherently a non-linear technique. The relations of many of the predictors, including El Nino itself with atmospheric circulation and rainfall is non linear (Drosdowsky and Williams 1991).
The forecast skill has been determined both for the rainfall principal components and on each individual grid point using the LEPS skill score. The properties of this measure are discussed in detail by Potts et al (1996).
A number of major problems with the forecast system, and any extension to other predictors and predictands, remain to be fully resolved.
1. The hindcast skill of the full SST scheme is reduced slightly in some regions and seasons, compared to the SOI or SST1 alone, by the inclusion of a number of additional predictors. This is particularly noticeable in the February-March-April season where there is now almost no skill anywhere compared with a few areas of weak to moderate skill using the SOI (or SST1 and SST2) alone, although the skill is generally improved through the "predictability barrier" from March-April-May onwards. This problem of decreased skill with increasing number of predictors, may be worsened if even more predictors other than the SST fields, eg indices of the atmospheric circulation such as the Quasi-Biennial Oscillation (QBO) are considered.
2. Each of the twelve overlapping seasons have been treated independently, and this has resulted, in some areas, in a lack of continuity in the selected predictors from season to season. In some circumstances this may cause major shifts in the forecast probabilities from season to season.
Both these problems, and the areas of large negative skill, are due to the large pool of potential predictors that the cross-model validation / selection procedure has to choose from. A partial solution to the first problem may be to pre-screen the field of potential predictors. This would probably need to be done separately for each season, as correlations between the SST and rainfall principal components suggest that all the SST principal components may be useful predictors in at least one season. This, however, may exacerbate the lack of continuity problem.
Finally, the question of the physical significance of some of the selected rainfall - SST relationships has not been examined. The links between SST1 (or the SOI) and SST2 and Australian rainfall are well documented. The relationships between Tasman Sea SSTa (SST9) and Australian rainfall have been explored (Streten, 1981, 1983; Priestly, 1965; Priestly and Troup, 1966), as have those with North Australian / Indonesian SSTa (SST12; Nicholls, 1984; Allan and Pariwono, 1990). Mid-latitude atmosphere - ocean interaction is a complex field, with most research (eg Frankingoul, 1985; Wallace and Jiang, 1987) suggesting that SSTa in this region are driven by atmospheric circulation anomalies. The observed mid-latitude SSTa patterns (eg SST2, SST4, SST5, SST6, SST7 and SST8) may then not be the direct cause of Australian rainfall variations, but instead may be indicative of large scale atmospheric circulation anomalies which are the true cause of the rainfall anomalies. This type of connection was demonstrated by Drosdowsky (1993b) for SST2 and SST6 and southern Australian winter rainfall. The connection between the North Pacific patterns and Australian rainfall may involve wave trains extending into both hemispheres from regions of equatorial convection, with only the Northern Hemisphere SSTs responding to the atmospheric anomalies. Further research is needed to clarify these relationships. ;
The research in this report was supported in part by a Land and Water Resources Research and Development Corporation (LWRRDC) grant, project No BOM1. Assistance and advice on various aspects of this project was provided by present and former members of the BMRC Climate Group; Neville Nicholls, Neville Smith, Carsten Frederiksen, David Jones, Richard Kleeman, Ramesh Balgovind and Alex Kariko. Constructive reviews and comments on this report were provided by David Jones, Bob Seaman, Mary Voice, John Cramb, Harvey Stern, Grant Beard and Huqiang Zhang.
Allan, R.J., and J.I. Pariwono, 1990: Ocean-atmosphere interactions in low-latitude Australasia. Intn. J. Climatol., 10, 145-178. ;
Barnston, A.G., 1994: Linear statistical short-term climate predictive skill in the Northern hemisphere. J. Climate, 7, 1513-1564. ;
Barnston, A.G., and H.M. van den Dool, 1993: A degeneracy in cross-validated skill in regression-based forecasts. J. Climate, 6, 963-977. ;
Bretherton, C.S., C. Smith, and J.M. Wallace, 1992: An intercomparison of methods for finding coupled patterns in climate data. J. Climate, 5, 541-560. ;
Casey, T., 1995: Optimal linear combination of seasonal forecasts. Aust. Meteor. Mag., 43, 219-224.
Drosdowsky, W., 1993a: An analysis of Australian seasonal rainfall anomalies: 1950-1987. I: Spatial patterns. Intl. J. Climatol., 13, 1-30.
Drosdowsky, W., 1993b: An analysis of Australian seasonal rainfall anomalies: 1950-1987. II: Temporal variability and teleconnection patterns. Intl. J. Climatol., 13, 111-149.
Drosdowsky, W., 1993c: Potential predictability of winter rainfall over southern and eastern Australia using Indian Ocean sea surface temperature anomalies. Aust. Meteor. Mag., 43, 1-6. ;
Drosdowsky, W., 1993d: Atlas of Australian monthly rainfall: 1910-1991. Unpublished BMRC Climate Report.
Drosdowsky, W., 1994: Analogue (non-linear) forecasts of the Southern Oscillation Index time series. Wea. and Forec., 9, 78-84.
Drosdowsky, W., and M. Williams, 1991: The Southern oscillation in the Australian region. Part I: Anomalies at the extremes of the oscillation. J. Climate, 4, 619-638.
Elsner, J.B., and C.P. Schmertmann, 1994: Assessing forecast skill through cross validation. Wea. and Forec., 9, 619-624.
Frankignoul, C., 1985: Sea surface temperature anomalies, planetary waves, and air-sea feedbacks in the middle latitudes. Rev. Geophys, 23, 357-390.
Ganeshanandam, S., and W.J. Krzanowski, 1989: On selecting variables and assessing their performance in linear discriminant analysis. Austral. J. Statist., 31, 433-447.
Jones, D.A. and G. Weymouth, 1997. An Australian monthly rainfall dataset. Technical Report No. 70, Bureau of Meteorology. Melbourne. 19pp.
Jones, D.A. and G. Beard, 1998; Verification of Australian Monthly District Totals using High Resolution Grided Analyses. Aust Meteor Mag. (accepted)
Michaelsen, J., 1987: Cross-validation in statistical climate forecast models. J. Climate Appl. Meteor., 26, 1589-1600.
McBride, J.L. and Nicholls, N., 1983: Seasonal relationships between Australian rainfall and the Southern Oscillation. Mon. Weath. Rev., 111, 1998-2004.
North, G.R., T.L. Bell, R.F. Cahalan, and F.J. Moeng, 1982: Sampling errors in the estimation of empirical orthogonal functions. Mon. Wea. Rev., 110, 699- 706. ;
Nicholls, N., 1984: The Southern Oscillation and Indonesian Sea surface temperature. Mon. Wea. Rev., 112, 424-432.
Nicholls, N., 1989: Sea surface temperatures and Australian winter rainfall. J. Climate, 2, 965-973.
Parker, D.E., C.K.Folland, A.C. Bevan, M.N. Ward, M. Jackson, and K.Maskell, 1995: Marine surface data for analysis of climatic fluctuations on interannual-to-century time scales. Natural Climate Variability on Decadal-to-Century Time Scales, National Research Council, 241-252.
Potts, J.M., C.K. Folland, I.T. Jolliffe, and D. Sexton, 1996: Revised "LEPS" scores for assessing Climate Model simulations and long-range forecasts. J.Climate, 9, 34-53.
Priestly, C.H.B., 1964: Rainfall - sea surface temperature associations on the New South Wales coast. Aust. Meteor. Mag., 47, 15-25.
Priestly, C.H.B., and A.J. Troup, 1966: Droughts and wet periods and their association with sea surface temperature. Aust. J. Sci., 29, 56-57.
Reynolds, R.W., 1982: Sea surface (and subsurface) temperatures for September 1980 - September 1981. Proceedings 6th Annual Climate Diagnostics Workshop, Palisades, New York, 14-16 October 1981.
Reynolds, R.W., and T.M. Smith, 1994: Improved global sea surface temperature analyses using optimal interpolation. J. Climate, 7, 929-948.
Smith, N.R., 1995; The BMRC ocean thermal analysis system. Aust. Met. Mag. 44, 93-110.
Stone, M., 1974: Cross-validatory choice and assessment of statistical predictions. J.R. Statist. Soc., B36, 111-147.
Stone, R.C., 1992: SOI phase relationships with rainfall in Australia. Ph.D. Thesis, The University of Queensland, St, Lucia, Queensland.
Streten, N.A., 1981: Southern hemisphere sea surface temperature variability and apparent associations with Australian rainfall. J. Geophys. Res., 86, 485-497.
Streten, N.A., 1983: Extreme distributions of Australian annual rainfall in relation to sea surface temperature. J. Climatol., 3, 143-153. ;
Troup, A.J., 1967: Opposition of anomalies in upper tropospheric winds at Singapore and Canton. Aust. Meteor. Mag., 15, 32-37.
Wallace, J.M., and Q. Jiang, 1987: On the observed structure of the interannual variability of the atmosphere/ocean climate system. Atmospheric and Oceanic Variability, Royal Meteorological Society, H.Cattle, Ed. 17-43.
Ward, N.M., and C.K. Folland, 1991: Prediction of seasonal rainfall in the north Nordeste of Brazil using eigenvectors of sea surface temperature. Intl. J. Climatol., 11, 711-743.
Zhang, X.-G., and T.M. Casey, 1992: Long-term variations in the Southern Oscillation and relationships with Australian rainfall. Aust. Meteor. Mag., 40, 211-225. ;
APPENDIX A: CROSS-VALIDATION ALGORITHM FOR A SINGLE PRESCRIBED MODEL
Given a data set of N cases, {xi, yi}, i=1,N ;
-----< Start loop over all N cases
|
| Set aside one case xk, yk
| Estimate model parameters
| from remaining (N-1) data,
| ie y = akx + bk
| Hindcast missing case yk
|
----> End loop over all N cases ;
Evaluate skill over the set of N independent hindcasts. This skill refers to the class of model, not the coefficients ak and bk.
To develop the forecast model for use on independent data use all N data points to fit best model y = ax + b ;
APPENDIX B: CROSS-VALIDATION ALGORITHM FOR MODEL SELECTION FROM MULTIPLE PREDICTORS
Given a data set of N cases, {x1i, x2i, ...., yi}, i=1,N
Specify pool of potential predictors, and class of models to select from
------< Start loop over all N cases
| Set aside one case x1k, x2k, ..., yk
| Find best model using remaining (N-1) cases
|
| ------< For each model,
| |
| | ------< Loop over all (N-1) cases
| | |
| | | Set aside one case x1j, x2j, ..., yj
| | | Estimate model parameters from remaining (N-2) data,
| | | ie y = a1jx1 + a2jx2 + ... + bj
| | | Hindcast missing case yj
| | |
| | ------> End loop over all (N-1) cases
| |
| | Evaluate skill of model over set of (N-1) independent hindcasts.
| |
| -------> End Loop over all Models
|
| Select model with most skill.
| Develop forecast model for use on independent set aside case k
| using all (N-1) cases to find parameters of best model
| ie y = a1kx1 + a2kx2 + ... + bk
| Hindcast missing case yk
|
--------> End loop over all N cases
Evaluate skill over the set of N independent hindcasts. This skill refers to the class of models evaluated, not a particular model or the coefficients of that model.
To develop the forecast model for use on independent data repeat the inner loops over all models but now using all N data points to find best model and model parameters
y = a1x1 + a2x2 + ... + b
Figure 1. Scree plot of eigenvalue versus eigenvector number for the first 24 eigenvectors of the correlation matrix of the standardised monthly anomalies of the two grided Australian rainfall data sets. Also shown are the estimated error bars for the eigenvalues based on the criterion of North et al (1982). Light trace is for the box averaged data, and heavy line for the analysed data set. The percentages of variance accounted for by the first component in each data set (21.7% for box average vs 22.7% for analysed) are comparable; the difference in eigenvalue is related to the different number of gridpoints (553 vs 631) used in each analysis.
Figure 2. Spatial pattern of loadings and associated scores (time series) of the first nine grided Australian rainfall VARIMAX rotated principal components of the standardised month anomalies of the data set. Contour interval is 0.2, with zero contour heavy, negative contours dashed and areas above +0.2 and below -0.2 shaded.
Figure 3. Scree plot of eigenvalue versus eigenvector number for the first 50 eigenvectors of the correlation matrix of the standardised monthly anomalies of GISST data set. Also shown are the estimated error bars for the eigenvalues based on the criterion of North et al (1982).
Figure 4. Spatial pattern of loadings and associated scores (time series) of the first twelve VARIMAX rotated principal components of the standardised month anomalies of the GISST data set. Contour interval is 0.2, with zero contour heavy, negative contours dashed and areas above +0.2 and below -0.2 shaded. First 6;; Second 6
Figure 5. Cross-validation LEPS scores for each of the twelve overlapping three monthly seasonal rainfall hindcasts using the previous season SOI (zero months lag), as predictor for the period 1950-1993.
Figure 6. As for Figure 5, but using the seasonal SOI with one month lag as predictor.
Figure 7. As for figure 5, but with single month SST1, lagged by one month as predictor.
Figure 8. As for Figure 5 but using bi-monthly values of the Indian Ocean Index, lagged by one and three months as predictors.
Figure 9. As for figure 5, but with single monthly values of SST2, lagged by one and three months as predictor.
Figure 10. As for figure 5, but with single monthly values of SST2, lagged by one month as predictor.
Figure 11. As for figure 5 but using both of the first two SST principal components (SST1 and SST2) lagged by one month as predictors.
Figure 12. As for Figure 5 but using the best subjectively determined combination of the first two SST principal components (SST1 or SST2 or both) or climatology as predictors.
Figure 13. As for Figure 5, but using the best, and possibly different, combination of the first two SST principal components or climatology, as determined by an inner cross validation predictor selection, as predictors.
Figure 14. Percentage of variance, calculated as the sum of the squared principal component loadings, explained by the first nine principal components of Australian monthly rainfall for 1950-1995.
Figure 15. Predictors to be used in independent forecasts for each rainfall principal component and season, selected by the cross validation technique applied to the full 1950-1993 data set. Principal components marked with an asterisk (*) are lagged by three months. SST12 in Figure 4 is denoted by 10, as SST10 and SST11 are not included in the pool of potential predictors.
Figure 16. "In sample" or single cross-validation LEPS scores for seasonal rainfall hindcasts using the best combination of one or two SST principal components selected from the full pool of twenty potential predictors for the period 1950-1993.
Figure 17. Independent, "out of sample" double cross-validation LEPS scores for seasonal rainfall hindcasts using the best combination of one or two SST principal components selected from the full pool of twenty potential predictors for the period 1950-1993.
Figure 18. LEPS skill scores for 42 independent hindcasts/forecasts for seasonal rainfall from JFM 1994 to JJA 1997. (a) Spatial pattern, and (b) temporal pattern averaged over each Australian state and the entire continent.
[an error occurred while processing this directive]
