|
||||||||||||||||||||||||||||||||||||||||||
BMRC is now part of CAWCR: The Centre for Australian Weather and Climate Research.
For more information on The Centre please go to http://www.cawcr.gov.au
|
BMRC Research Report No 65
NEAR GLOBAL SEA SURFACE TEMPERATURE ANOMALIES AS PREDICTORS OF AUSTRALIAN SEASONAL RAINFALL
Wasyl Drosdowsky
and
Lynda E. Chambers
January 1998
Contents
4 Principal Component Analysis
(c) Full SST data set - multiple predictors.
Acknowledgments
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.
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.
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.
4. PRINCIPAL COMPONENT ANALYSIS
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).
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.
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.
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.
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".
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.
Acknowledgments
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.
| Bureau Home || BMRC Home || Search || Contact BMRC Webmaster |
| Experimental results described in these pages are from research systems developed in BMRC and are not part of the Bureau of Meteorology's operational products & services. |
Home | About Us | Learn about Meteorology | Contacts | Search | Help | Feedback Weather and Warnings | Climate | Hydrology | Numerical Prediction | About Services | Registered Users | SILO |
|
© Copyright Commonwealth of Australia 2008, Bureau of Meteorology (ABN 92 637 533 532) Please note the Copyright Notice and Disclaimer statements relating to the use of the information on this site and our site Privacy and Accessibility statements. Users of these web pages are deemed to have read and accepted the conditions described in the Copyright, Disclaimer, and Privacy statements. Please also note the Acknowledgement notice relating to the use of information on this site. No unsolicited commercial email. |
