Modelling current trends in Northern ... - Wiley Online Library

INTERNATIONAL JOURNAL OF CLIMATOLOGY Int. J. Climatol. 26: 867–884 (2006) Published online 26 January 2006 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/joc.1286

MODELLING CURRENT TRENDS IN NORTHERN HEMISPHERE TEMPERATURES TERENCE C. MILLS* Department of Economics, Loughborough University, Loughborough, LE11 3TU, UK Received 9 August 2005 Revised 12 October 2005 Accepted 17 October 2005

ABSTRACT Fitting a trend is of interest in many disciplines, but it is of particular importance in climatology, where estimating the current and recent trend in temperature is thought to provide a major indication of the presence of global warming. A range of ad hoc methods of trend fitting have been proposed, with little consensus as to the most appropriate techniques to use. The aim of this paper is to consider a range of trend extraction techniques, none of which require ‘padding’ out the series beyond the end of the available observations, and to use these to estimate the trend of annual mean Northern Hemisphere (NH) temperatures. A comparison of the trends estimated by these methods thus provides a robust indication of the likely range of current trend temperature increases and hence inform, in a timely quantitative fashion, arguments based on global temperature data concerning the nature and extent of global warming and climate change. For the complete sample 1856–2003, the trend is characterised as having long waves about an underlying increasing level. Since around 1970, all techniques display a pronounced warming trend. However, they also provide a range of trend functions so that extrapolation far into the future would be a hazardous exercise. Copyright  2006 Royal Meteorological Society. KEY WORDS:

Northern Hemisphere temperature; trends; trend extraction; filters; smooth transition; non-parametric modelling

1. INTRODUCTION Fitting a trend is of interest in many disciplines, such as economics and finance (see Mills, 2003, for a survey of techniques in these areas). It is of particular importance, however, in climatology, where estimating the current and recent trend in temperature is thought to provide a major indication of the presence of global warming. A range of ad hoc methods of trend fitting have been proposed, with little consensus as to the most appropriate techniques to use. Among recent contributions are Benestad (2003), Jones and Lister (2004), Moberg and Jones (2005) and Schmidlii and Frei (2005), who all fit deterministic time trends, either linear, polynomial or logistic, Piˇsoft et al. (2004), who use wavelet transforms, and Brázdil et al. (2005), who employ low-pass filters. A pertinent example of this lack of consensus is the current debate on how best to use filters to compute current trend estimates for Northern Hemisphere (NH) air and sea-surface temperatures, a crucial input into debates concerning global warming: see, e.g. Houghton et al. (2001) and, for a ‘sceptical’ critique, Essex and McKitrick (2002). Rather than arguing over the form of trend filter that is used to compute such estimates, and hence the implicit underlying trend model, this debate has focused on the best way of extending the series – of ‘padding’ it out – so that a two-sided filter can be used right up to the end of the sample: see, amongst others, Mann and Jones (2003), Mann (2004) and Soon et al. (2004). The aim of this paper is to consider a range of trend extraction techniques, none of which require ‘padding’ out the series beyond the end of the available observations, and to use these to estimate the trend of annual * Correspondence to: Terence C. Mills, Department of Economics, Loughborough University, Loughborough, LE11 3TU, UK; e-mail: [email protected]

Copyright  2006 Royal Meteorological Society

868

T. C. MILLS

mean NH temperatures. A comparison of the trends estimated by these methods will thus provide a robust indication of the likely range of current trend temperature increases and hence inform, in a timely quantitative fashion, arguments based on global temperature data concerning the nature and extent of global warming and climate change. Note that, in line with the debate mentioned earlier, the trend extrapolation techniques considered here are all univariate, and thus use only the history of the NH temperature series and not the information provided by possible ‘covariates’. The methods include two-sided filters whose weights are adjusted near the end of the sample to enable them to be used for calculating current trend estimates without the need to invent future values.

2. TECHNIQUES FOR TREND EXTRACTION The definition of ‘trend’ has long been a vexed question for statisticians. An apposite view is that of Kendall et al. (1983), §45.12, who muse that ‘generally, one thinks of [the trend] as a smooth broad motion of the system over a long term of years, but ‘long’ in this connexion is a relative term, and what is long for one purpose may be short for another. For example, if we were examining rainfall records over a hundred years, a slow rise from the beginning of the period to the end would be regarded as a trend; but if we possessed records for two thousand years . . . the rise over a particular century might appear as a slow oscillatory movement, so that any inference from the ‘trend’ in a particular century to the effect that the weather was likely to continue becoming wetter might be quite false.’ Furthermore, they go on to say that, ‘if we were attempting to study climatic changes over the face of the earth for geological periods of time we should accept the continuance of the trend with the greatest reserve or, more probably, reject it on collateral grounds.’ Harvey (2002, page 2243) is more precise, viewing the trend as ‘that part of the series which when extrapolated gives the clearest indication of future long-term movements in the series.’ These twin ideas of smoothness and prediction inform our discussion and choice of trend extraction techniques. Equally as important, though, is the ability of the technique to estimate a trend for the complete sample, particularly those dates close to the end of the available record, without the necessity of extrapolating the series into the future in some fashion. Six methods are thus investigated, which it is convenient to group into four categories: ARIMA decomposition, parametric and non-parametric trends, and low-pass filters. All, however, rely on the basic decomposition of an observed series, Zt , into unobserved trend and, for want of a better term, noise components, St and Nt , respectively: Zt = St + Nt

t = 1, 2, . . . , T

(1)

To focus ideas on the problem at hand, Figure 1 plots annual mean NH temperatures from 1856 to 2003. This series is the focus of the debate alluded to in the ‘Introduction’ and is obtained from the Climatic Research Unit (CRU) of the University of East Anglia, UK, at: http://www.cru.uea.ac.uk/ftpdata/tavenh2v.dat; see Jones et al. (1999) and Jones and Moberg (2003) for detailed discussion. This plot reveals the presence of changing trend behaviour throughout the sample period, with fairly homogeneous fluctuations about trend. The behaviour of the series over the last three decades – a distinct upward movement after a similar period of relative stability – is the obvious focus of attention in debates concerning global warming, and the estimation of this recent trend is of primary concern here. 2.1. ARIMA decomposition A popular method of estimating the current and recent trend in a time series is to compute the Beveridge and Nelson (BN;1981) trend component from the ARIMA representation of Zt . Suppose Zt is integrated of order one (denoted I (1)) and that Zt = (1 − B)Zt = Zt − Zt−1 is generated by an ARMA (p,q) process, i.e. Zt = β +

(1 + θ1 B + · · · + θq B q ) θ (B) at = β + at (1 + φ1 B + · · · + φp B p ) φ(B)


(2) Int. J. Climatol. 26: 867–884 (2006)

869

MODELLING CURRENT TRENDS IN NH TEMPERATURES

Temperature (degrees C) relative to 1961-1990 mean

.6 .4 .2 .0 -.2 -.4 -.6 1875

1900

1925

1950

1975

2000

Year

Figure 1. Annual mean Northern Hemisphere temperatures: 1856–2003

The trend component St can then be estimated recursively as (see Appendix A1 for details) St = −θ1 St−1 − · · · − θq St−q + (θ (1)/φ(1))(Zt + φ1 Zt−1 + · · · + φp Zt−p ), t > max(p, q)

(3)

St = Zt ,

(4)

t ≤ max(p, q)

The BN trend is thus estimated as a weighted average of the past trend and observed values, with the noise component being obtained by residual. It is thus a one-sided filter and hence does not require future values of Zt to estimate the current trend St . 2.2. Parametric techniques The following techniques, structural time series models and smooth transitions, are labelled parametric as they both propose formal models for the trend component, St , whose parameters may be estimated from the observed data. They differ in that the former allows for a trend component that may be stochastic, while the latter forces the trend to follow a deterministic function of time, both being assumed to be independent of the noise Nt . Structural models have been used to extract unobserved components from the North Atlantic Oscillation (Mills, 2004) and from the England and Wales rainfall data (Mills, 2005). Applications of structural models to estimate trends in temperatures and environmental series may be found in Stern and Kaufmann (2000) and Visser (2004) respectively. Smooth transition models have been fitted to Northern Hemisphere, Southern Hemisphere and global temperature data by Harvey and Mills (2001); see also Harvey and Mills (2002). 2.2.1. Structural time series models. Figure 1 shows that the NH temperature series may be characterised as going through a sequence of trend regimes, moving relatively smoothly from one regime to the next. It may thus be said to exhibit local trend behaviour. The basic structural model (BSM; see e.g. Harvey (1989) for extended discussion and development) allows such behaviour to be modelled by assuming that the trend follows a random walk with a stochastic drift (slope), which also follows a random walk: St = St−1 + βt−1 + v1t

(5)

βt = βt−1 + v2t

(6)

The level and slope innovations, v1t and v2t , are assumed to be mutually uncorrelated zero mean white noises. Several special cases of the BSM often occur in practice. If both the level and slope innovation variances, σ12 and σ22 , are zero, then βt = βt−1 = β and St will be the deterministic linear trend St = µ + βt Copyright  2006 Royal Meteorological Society

(7) Int. J. Climatol. 26: 867–884 (2006)

870

T. C. MILLS

where µ = S0 can be interpreted as the ‘initial condition’ for the trend. If just σ22 = 0, the slope is fixed and the trend reduces to a random walk with constant drift: St = St−1 + β + v1t

(8)

If, however, just σ12 = 0, then the trend becomes 2 St = v2,t−1

(9)

which is an integrated random walk and is often referred to as the ‘smooth trend’ model since fitted models of this type often produce a trend that is relatively smooth. In general, the noise component may be defined as the sum of a stationary autoregressive component, a stochastic cyclical component and a white noise ‘irregular’ component. As these are not our focus of interest here, their exact structure need not concern us, but Koopman et al. (2000) may be consulted for a general treatment of structural models of this type. The unknown parameters of the model and the various components may be estimated by casting the model into state space form and employing the Kalman filter. Maximum likelihood estimates of the parameters, typically the innovation variances, may then be obtained using the EM algorithm, which is a recursive method of obtaining maximum likelihood estimates of the unknown elements in the state and error system matrices of the state space form. Koopman (1993) provides a detailed development of the EM algorithm for the present set-up, and estimates of the components can be calculated by using a smoothing procedure. This first estimates the components recursively for t = 1, 2, . . . , T , thus obtaining the ‘filtered’ components, and then runs a second ‘backward’ recursion over t = T , T − 1, . . . , 1 to obtain the ‘smoothed’ estimates. Full details of the estimation procedure, incorporated in the software package STAMP, may be found in Koopman et al. (2000). 2.2.2. Smooth transition trends. An alternative approach to modelling the smoothly changing local trend behaviour that is observed in Figure 1 is to use a deterministic nonlinear trend. The logistic smooth transition function has a long tradition in the statistical modelling of changing regimes, being introduced by Bacon and Watts (1971) and extended to time series models by Teräsvirta (1994) and Leybourne et al. (1998). Harvey and Mills (2001, 2002) consider a general smooth transition trend model that allows for the possibility of three different regimes, taking the form St = α0 + β0 t + (α1 + β1 t)θ1t (γ1 , τ1 ) + (α2 + β2 t)θ2t (γ2 , τ2 )

(10)

where the θit (γi , τi ) are logistic smooth transition functions controlling the transition between regimes, defined by θit (γi , τi ) = [1 + exp{−γi (t − τi T )}]−1

i = 1, 2

(11)

The midpoints of the two transitions are given by τ1 T , and τ2 T for sample size T it is easily seen that for positive values of γi , θi,−∞ (γi , τi ) = 0. θi,∞ (γi , τi ) = 1 and θi,τi T (γi , τi ) = 0.5. The transition speeds are given by γ1 and γ2 . Small values of these parameters mean that the transition functions take a long time to traverse the interval [0,1] – in the limiting case when γi = 0, θit (γi , τi ) = 0.5 for all t and there is no transition. For large values, the transition functions traverse the interval [0,1] very quickly and as γi → ∞, θit (γi , τi ) changes from 0 to 1 instantaneously at time τi T . The transition paths are symmetric about their respective midpoints, although this could be relaxed as in Sollis et al. (1999). The three regimes admitted by the trend function in Equation (10) are given by α0 + β0 t, (α0 + α1 ) + (β0 + β1 )t, and (α0 + α1 + α2 ) + (β0 + β1 + β2 )t, with the two transition paths linking the three regimes allowed to differ in speed. Such a specification, if found to be acceptable, is capable of modelling the stylised facts of NH temperature data seen in Figure 1. These show generally stable temperatures up to the turn of the century, then a short decline followed by an increasing trend up to 1940, a period of slightly declining temperature until 1965, and finally a third warming period after this date. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)


871

Estimation of the smooth transition model can be viewed as an exercise in nonlinear least squares regression of the equation Zt = α0 + β0 t + (α1 + β1 t)θ1t (γ1 , τ1 ) + (α2 + β2 t)θ2t (γ2 , τ2 ) + Nt

(12)

where the noise Nt may, in general, be allowed to follow a stationary ARMA process, independent of St . Hypothesis tests may then be performed on the estimated coefficients if simpler trend functions suggest themselves, so that restricted smooth transition trends may end up being fitted. For example, insignificant α2 and β2 estimates would imply that only a single transition (two regime) model was required. Although more than three regimes could be considered, successfully identifying the regimes would be extremely difficult with a series of this length. 2.3. Non-parametric regression Suppose we wish to fit a very general trend function St = f (t) to Zt , i.e. we consider the model Zt = f (t) + Nt

(13)

A popular choice for f (t) is a local polynomial so that at each point in time τ the parameters of the local pth order polynomial are chosen to solve the weighted least squares problem Qp =

T

(Zt − δ0 − δ1 (τ − t) − · · · − δp (τ − t)p )2 K(τ − t; b)

(14)

t=1

where K(τ − t; b) is a symmetric positive function, often called a kernel, that decreases as |τ − t| increases. Here, b > 0 is a smoothing parameter, usually referred to as the bandwidth. The larger b is, the ‘greater’ the degree of averaging, and the smoother will be the estimated function fˆb,p (τ ). As is clear by this notation, the parameters p and b must be chosen by the analyst. With regard to the choice of polynomial order, theoretical research has shown that even values of p lead to greater asymptotic bias in fˆb,p (τ ) near the ends of the sample than do odd values. This asymptotic bias is a consequence of having to truncate the kernel weights at the end of the sample. Since accuracy at the end of the sample is important here, this suggests that settings of p = 1 or 3 should be preferred over p = 0 or 2. The bandwidth can be chosen through cross validation. This uses the idea of ‘leave-one-out’ prediction and minimises the criterion C(b) =

T 1 (τ ) (Zt − fˆb,p (t))2 T t=1

(15)

(τ ) (t) is the estimate of fˆb,p (t) based on the sample Z1 , Z2 , . . . , Zτ −1 , Zτ +1 , . . . , ZT , i.e. with where fˆb,p observation τ removed. Local polynomial fitting by kernel regression chooses K(τ − t; b) to be a density function (more precisely a function that integrates to unity). Although a variety of kernel functions are available, a popular choice is the Gaussian kernel, which, for a given choice of bandwidth b, weights the local observations according to

1 K(τ − t; b) = √ exp(−(τ − t)2 /2b2 ) b 2π

(16)

Thus, b plays the role of ‘standard deviation’ and hence the larger b is, the higher are the weights given to observations further from τ , and hence smoother the fitted trend function will be. Note that setting K(τ − t; b) = 1 yields the local polynomial trend analysed by Kendall (1973). Non-parametric regression is discussed in, for example, Härdle (1990), Simonoff (1996), Mills (2003) and Ruppert et al. (2003) and has been used by Hall and Tajvidi (2000) and Harvey and Mills (2003) to model various aspects of temperature trends. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)

872

T. C. MILLS

2.4. Hodrick–Prescott and Butterworth low-pass filtering techniques A popular trend extraction technique is to estimate the signal in Equation (1) by the n = 2r + 1 moving average filter r

Sˆt(r) =

wj Zt+j

(17)

j =−r

The weights in the moving average can be chosen in various ways: a ‘low-pass’ filter sets the weights of the moving average as wL,0 =

λc , π

wL,j =

sin(λc j ) πj

j = 1, . . . , r

(18)

Here λc is the ‘cut-off frequency’ such that when analysed in the frequency domain, the filter only passes frequencies in the interval −λc < λ < λc . Filters of this type were first used by Craddock (1957) to model temperature trends. However, choosing r is not simple, as too small a value leads to approximation problems known as leakage, compression, and exacerbation, while setting r too large leads to more ‘lost’ observations at the beginning and end of the sample, exactly the problem we wish to avoid here. Because of these difficulties, various other weighting mechanisms have been proposed, which attempt to provide good approximations to the ideal low-pass filter and, moreover, deal with the finite sample problem inherent in Equation (17). An approximate low-pass filter commonly used in economics is that associated with Hodrick and Prescott (1997). This filter, generally known as the HP (trend) filter, was originally developed as a solution to the problem of minimising the variation in the noise component subject to a smoothness condition that penalises acceleration in the trend. This is an approach that has a long history of use in, for example, actuarial science (Whittaker, 1923). Appendix A8 shows that the HP trend can be written as St =

1 Zt 1 + ξ(1 − B)2 (1 − B −1 )2

(19)

where ξ is a parameter that controls the smoothing of the trend. For an infinite sample, the HP filter thus estimates the trend as SˆtHP =

∞

wj Zt+j = wHP (B)Zt

(20)

j =−∞

where wHP (B) =

1 1 + ξ(1 − B)2 (1 − B −1 )2

(21)

Appendix A8 sets out how the filter is amended for a finite sample. An alternative trend filter estimate is wB (B) =

1+ξ

1 (1 − B)(1 − B −1 ) (1 + B)(1 + B −1 )

f

(22)

The filter wB (B) is known as the Butterworth (square-wave) filter and it has a long history of use in engineering (Pollock, 1999, 2000 and the references therein). Although wB (B) is a two-sided infinite filter, Pollock (2000) presents an algorithm that enables the trend filter to be calculated in finite samples. The smoothing parameter ξ in the HP and Butterworth filters may be chosen to correspond to a frequency cut-off of λc : Appendix A8 shows how. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)

873


2.5. A comparison of the trend models In summarising the similarities and differences between the various trend models, we make the following comments. The smooth transition and non-parametric models are regression-based approaches and thus require only observed data. The trend and noise components in these models are independent by construction. All parameters in the smooth transition are estimated, while non-parametric regression requires choosing the order of the local polynomial, the kernel, and the bandwidth, which control the smoothness of the fitted trend. The components in the BSM are independent by assumption and are calculated by an implied one-sided filter of current and past observed values whose weights depend on the estimated innovation variances. In contrast, the BN components are perfectly correlated, but again they are estimated using a one-sided filter of current and past observed values, the weights of which are functions of the estimated ARIMA parameters fitted to the data. The low-pass filters use asymmetric-moving averages that adjust the weights as the end of the series is approached. Each requires choosing various parameters that will control the smoothness of the fitted trends. For only the parametric BN, BSM, and smooth transition methods can standard errors of the trends be calculated, as these have formally specified trend components. Non-parametric regression and low-pass filtering attempt to approximate an unknown trend function and thus measures of precision are not available.

3. ESTIMATING TRENDS FOR NH TEMPERATURES 3.1. BN decomposition Figure 2 shows the sample autocorrelation function (SACF) for the levels (Zt ) of the annual mean NH temperature series for lags up to k = 36. The slow decline of the SACF suggests that Zt is non-stationary, and this is confirmed by a battery of unit root tests. For example, the Augmented Dickey–Fuller regression, with lag augmentation selected using Schwarz’s Information Criterion, is Zt = 0.011 − 0.010 Zt−1 + (0.012) (0.051)

3

δî Zt−i + et

(23)

i=1

Standard errors are shown in parentheses so that the Dickey–Fuller (1979) statistic (the t-statistic on the coefficient of Zt−1 ) for testing the null of a unit root in Zt is just −0.20, which has a probability value of 0.93 and therefore cannot reject this null. A similar test on Zt rejects the hypothesis that there is a unit root in the first differences so that Zt is indeed I (1) and can be modelled by an ARIMA (p,1,q) process. Figure 3 shows the SACF and sample partial ACF (SPACF) for Zt : the effective ‘cut-off’ of the SACF at lag 4 and the slow decline of the SPACF suggest that a moving average of order 4 is likely to provide an acceptable fit to the series. Estimation leads to the (restricted, since the constant and at−2 are able to be excluded as they 1.0 0.8

rk

0.6 0.4 0.2 0.0 5

10

15

20

25

30

35

k

Figure 2. SACF for annual mean Northern Hemisphere temperatures, 1856–2003 (Zt ). rk denotes the sample autocorrelation of Zt at lag k. The dashed horizontal line indicates an approximate two standard error bound (0.164) Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)

874

T. C. MILLS (a) 1.0

rk

0.5 0.0 -0.5 -1.0 5

10

15

20

25

30

35

20

25

30

35

k (b) 1.0

fkk

0.5 0.0 -0.5 -1.0 5

10

15 k

Figure 3. SACF (a) and SPACF (b) of the first differences of mean Northern Hemisphere temperatures: 1857–2003 (Zt ). rk denotes the sample autocorrelation of Zt at lag k. The sample partial autocorrelation φkk is the estimate of the coefficient on Zt−k in an AR(k) model fitted to Zt . Dashed horizontal lines indicate an approximate two standard error bound (0.164)

are statistically insignificant) ARIMA (0,1,4) model Zt = at − 0.604 at−1 − 0.309 at−3 + 0.265 at−4 (0.069) (0.082) (0.082)

σˆ a = 0.129

(24)

Diagnostic tests for residual autocorrelation, heteroskedasticity, and non-normality all prove insignificant so that this model provides a valid base for calculating the BN decomposition. Since φ(B) = 1 and θ (B) = 1 − 0.604B − 0.309B 3 + 0.265B 4 , φ(1) = 1, θ (1) = 0.352 and thus St = 0.604St−1 + 0.309St−3 − 0.265St−4 + 0.352Zt ≈ 0.65(0.9St−1 + 0.5St−3 − 0.4St−4 ) + 0.35Zt St = Zt

t = 5, . . . , T

(25)

t = 1, . . . , 4

(26)

Hence, the current value of the BN trend component is calculated as a weighted average (with weights (0.35,0.65)) of the current value Zt and a weighted average of the previous four trend values, these having weights (0.9,0,0.5,−0.4). It is thus a generalisation of the recursive form of simple exponential smoothing (Newbold and Vougas, 1996). It might be argued that the pattern of the SACF of Zt shown in Figure 2 could be better characterised as a hyperbolic decline, rather than the linear decline associated with I (1) non-stationarity. This would suggest that Zt requires fractional differencing to induce stationarity so that an ARFIMA (p, d, q) model, with a non-integer d, may be appropriate. In fact, allowing d to be non-integer obtained 1.03 Zt = at − 0.628 at−1 − 0.307 at−3 + 0.274 at−4 (0.142) (0.083) (0.094)

σˆ a = 0.128

(27)

As the standard error attached to dˆ = 1.03 is 0.12, d is insignificantly different from unity, and thus we do not consider the ARFIMA model as the basis for constructing trends. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)


875

The BN trend is shown superimposed upon the NH temperature series in Figure 4. Although the trend is considerably smoother than the observed data, it still displays noticeable short-term variation, as it is a one-sided filter throughout the sample period. 3.2. Parametric trends The best fitting BSM was found to contain the smooth trend 2 St = νt

σ1 = 0

σˆ 2 = 0.00361

(28)

or, in structural form, St = βt βt = βt−1 + νt

(29)

The noise component Nt was found to be satisfactorily modelled as the sum of a white noise irregular component and a cyclical component having a stochastic cycle with an average period of seven years, which we regard as too short a cycle to be part of the trend. As the variance of Nt is σˆ N2 = 0.10642 , the ‘signal-tonoise’ variance ratio is σˆ 22 /σˆ N2 = 0.0012, so that the trend is indeed smooth, as can be seen from Figure 5. The smooth trend has local minima at 1910 and 1971 and ‘internal’ local maxima at 1873 and 1946. The smooth transition trend was estimated as St = − 0.274 + − 1.133 + 0.0179 t θ1t (0.56, 0.34) (0.018) (0.325) (0.0053) 0.0147 2.786 t θ2t (0.11, 0.73) + + − (30) (1.056) (0.0066) Estimation was carried out by ordinary least squares as the residuals from the regression (effectively the estimates of Nt ) were found to be indistinguishable from Gaussian white noise. The parameter β0 was estimated to be insignificantly different from zero and thus the trend term in the first regime was set to zero. The three trend regimes are therefore estimated to be −0.274, −1.407 + 0.0179t, and −4.193 + 0.0326t, with transition midpoints estimated at 1906 and 1965, although the local minima are slightly later at 1909 and


.6 .4 .2 .0 -.2 -.4

BN

-.6 1860 1870 1880 1890 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 Year

Figure 4. Beveridge–Nelson trend Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)

876

T. C. MILLS


.6 .4 .2 .0 -.2 -.4

BSM Smooth transition

-.6 1875

1900

1925 Year

1950

1975

2000

Figure 5. Basic structural model and smooth transition trends

1970. The internal local maximum is again at 1946. The smooth transition trend is also plotted in Figure 5 and differs from the BSM trend primarily by forcing St to be flat during the 1800s. 3.3. Non-parametric trend A local cubic (p = 3) trend is shown in Figure 6, having been computed by setting the bandwidth at b = 12, this being determined by cross validation. Setting p at a lower value produced a less smooth trend. The trend has local minima at 1908 and 1970 and local maxima at 1871 and 1948. 3.4. Low-pass filters The spectral density of the NH temperature series is shown in Figure 7. From this, a low-pass frequency cut-off of λc = π/16 ≈ 0.2 was chosen, which, from (A32) and A(33), sets the smoothing parameters as ξH P = 3020, and with f = 4, ξB = 9.74 × 106 . This choice of f is one of those recommended by


.6 .4 .2 .0 -.2 -.4

Local cubic

-.6 1860 1870 1880 1890 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 Year

Figure 6. Local cubic polynomial trend Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)


877

.22 .20 .18 .16

Spectral density

.14 .12 .10 .08 .06 .04 .02 .00 π /16

π /2

π

Frequency

Figure 7. Spectral density of mean annual NH temperature; λc = π/16

Pollock (2000). The setting of λc , which represents a compromise between capturing as many low frequency components as possible while retaining smoothness, is rather higher than the value of 0.025 chosen arbitrarily by the protagonists in the global warming debate for their (rather different) low-pass filters (e.g. Mann, 2004, and Soon et al., 2004). The HP and Butterworth trends are shown in Figure 8 and, apart from some differences at the start of the sample, are close to each other and, indeed, to the non-parametric trend shown in Figure 6. Both low-pass trends date the late local minimum at 1971 and the early one at 1910 for the HP and at 1911 for the Butterworth trend. The HP trend dates the local maxima at 1873 and 1948, but the Butterworth trend puts these as 2 years earlier. 3.5. Modelling recent trends in NH temperatures Although there are some differences, most notably that the BN trend is much less smooth than the others, the estimated trends shown in Figures 4–6 and 8 exhibit marked ‘long wave’ patterns, with peaks around 1873 and 1948 and troughs around 1909 and 1970. Since this last trough, the temperature trend has been steadily rising and it is this period that has naturally attracted attention in debates concerning global warming. The rate of trend increase in the period since 1970 has been twice as fast as that of the earlier warming period between 1910 and 1948. This period saw average annual trend increases ranging from 0.009 to 0.012 ° C across the six estimated trends, whereas the most recent period has had annual trend increases from 0.018 to 0.021 ° C. We now focus on the recent period from 1993 to 2003. Excluding the BN decomposition, trend temperature levels range from 0.305 to 0.313 ° C in 1993, whereas in 2003 the levels range from 0.575 to 0.613 ° C, reflecting the manner in which the alternative trend models respond to the most recent observations at the end of the sample. This is clearly observed in Figure 9, which shows the annual trend increases over this period. Trend increases in 2003 range from 0.0275 to 0.0315 ° C, although it is also seen from this figure that the rate of trend increase appears to be declining. The BN trend is omitted from Figure 9 as the variability of its values disrupts the scaling: over this recent period BN trend increases had a maximum value of 0.118 in 1998 and Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)

878

T. C. MILLS


0.6 0.4 0.2 0.0 -0.2 -0.4

HP Butterworth

-0.6 1875

1900

1925

1950

1975

2000

Year

Figure 8. HP and Butterworth trends

0.032 Smooth transition Local cubic

0.031

Temperature change

0.030 BSM

0.029

Butterworth

0.028

HP 0.027

0.026

0.025 1994

1996

1998

2000

2002

Year

Figure 9. Estimated trend increases using the various trend estimation techniques; 1993–2003

a minimum of −0.041 in 1996. Its trend level is also always below all the other trends in this period. This is because the omission of the (statistically insignificant) constant from the ARIMA representation implies that the BN trend evolves as a driftless random walk, so that it has difficulty in adapting to the change in the trend seen in the observed series since 1970. It is a little difficult to compare these latest trend estimates with those calculated by the protagonists in the debate alluded to in the ‘Introduction’, as a range of end years from 2000 have been used in the various papers. The closest comparison we can make is with Mann (2004), whose preferred trend estimate for 2003 is 0.610 ° C. This is obtained using a ‘Butterworth-type’ low-pass filter with λc = 0.025 and pads the series out into the future ‘with values within one filter width of the boundary reflected vertically (i.e. about the ‘y’ axis) Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)


879

relative to the final value. This tends to impose a point of inflection at the boundary, and leads the smooth towards the boundary with constant slope’ (Mann, 2004, page 1). This estimate is at the top end of the range of trend estimates provided here. By contrast, two other methods of data padding used, but not favoured, by Mann produce trend estimates of approximately 0.4 and 0.2 ° C for 2003, suggesting that trends estimated by these methods are essentially arbitrary.

4. DISCUSSION AND SUMMARY A range of trend extraction techniques has been examined in the context of estimating the trend in annual NH temperatures. Their common property is that they do not require extrapolating data into the future to be able to estimate trend values at the end of the sample. These techniques are by no means exhaustive of those available that have this property. For example, Bell and Martin (2004) provide a method for estimating the trend when the underlying unobserved components follow ARIMA processes, while Gray and Thomson (2002) and Quenneville et al. (2003) extend the Henderson–Musgrave-type trend filters used in seasonal adjustment procedures such as X-12-ARIMA (Findlay et al., 1998) to enable them to be used right up to the end of the sample. The technique that stands apart from the others is the BN trend, which is calculated from the ARIMA representation of the series and is not based on any smoothing of the data. This feature alone would set the BN trend apart, but its singularity is exacerbated by it being, in effect, a driftless random walk, so that it captures poorly the increase in temperatures that have occurred since around 1970. The low-pass filters provide pleasingly smooth representations of the trend throughout the sample but also the two flattest trends and the smallest trend increases over the last ten years. Adjusting the smoothness parameters (equivalently, changing the frequency cut-off in Figure 7) led to trends that were either far too smooth, thus missing the ‘long waves’ in temperature, or to trends that were too ‘wiggly’. The parametric trends have an advantage in that estimates of the precision with which the trend and its ‘slope’ are estimated can be calculated. The BSM smoothing algorithm estimates the 2003 trend to be 0.595 ° C with a standard error of 0.057, with the associated slope, βˆ2003 = 0.029, having a standard error of 0.010 (see Koopman et al., 2000, for details of the algorithm used to compute such standard errors). For the smooth transition model, the 2003 trend estimate is 0.613 with a standard error of 0.115. The slope of the trend in the final regime is estimated to be 0.033 and has a standard error of 0.008 attached to it. Interestingly, since the 2003 slope for this model is estimated as 0.031, the series is still in transition from the second to the third regime. Nonetheless, these estimates of precision are such that standard confidence intervals for the trend slope would contain all six estimated slopes. (The standard error of the BN slope is 0.097, reflecting the much greater variability of this trend estimate.) In summary, then, although the techniques investigated here display certain idiosyncrasies, they share the common feature of a ‘long wave’ in trend NH temperatures with a pronounced warming trend since 1970. The range of trend functions provided by these techniques should, however, temper the enthusiasm of anyone tempted to use them for extrapolating trend temperatures far into the future.

APPENDIX A.1. The Beveridge–Nelson decomposition Suppose Zt is integrated of order one (denoted I (1)) and that Zt = (1 − B)Zt = Zt − Zt−1 has the Wold decomposition: Zt = β +

∞

ψj at−j = β + ψ(B)at

ψ0 = 1

(A1)

j =0


Int. J. Climatol. 26: 867–884 (2006)

880

T. C. MILLS

Here, ψ(B) = 1 + ψ1 B + ψ2 B 2 + · · · is an infinite order polynomial in the lag operator B, defined as B at ≡ at−j , and the innovation at is zero mean white noise with variance σa2 , so that E(at at−j ) = 0 for all 2 2 j = 0 and E(at ) = σa . Since ψ(1) = ψj is a constant, we may write, in general, j

ψ(B) = ψ(1) + C(B)

(A2)

so that C(B) = ψ(B) − ψ(1) = −ψ1 (1 − B) − ψ2 (1 − B 2 ) − ψ3 (1 − B 3 ) − · · · = (1 − B)(−ψ1 − ψ2 (1 + B) − ψ3 (1 + B + B 2 ) − · · ·)         ∞ ∞ ∞ = (1 − B) −  ψj  −  ψj  B 2 −  ψj  B 3 − . . . j =1

j =2

j =3

˜ = ψ(B)

(A3)

Thus, ˜ ψ(B) = ψ(1) + ψ(B)

(A4)

so that ˜ Zt = β + ψ(1)at + ψ(B)a t

(A5)

and, in terms of the decomposition Zt = St + Nt , St = β + ψ(1)at

(A6)

˜ Nt = ψ(B)a t

(A7)

and

Thus the Beveridge–Nelson (BN) trend is a random walk with drift β and innovation ψ(1)at , which will therefore be proportional to the innovation to Zt , and will have a variance of (ψ(1))2 σa2 , which may be larger or smaller than σa2 depending on the signs and patterns of the ψj s. The noise component is clearly stationary, but since it is driven by the same innovation as the trend component, St and Nt must be perfectly correlated. If the Wold decomposition Equation (A1) is approximated by an ARMA (p,q) process, then Zt = β +

(1 + θ1 B + · · · + θq B q ) θ (B) at p at = β + (1 + φ1 B + · · · + φp B ) φ(B)

(A8)

St = β +

(1 + θ1 + · · · + θq ) θ (1) at = β + at (1 + φ1 + · · · + φp ) φ(1)

(A9)

and


Int. J. Climatol. 26: 867–884 (2006)


881

Now, Equation (A8) can be written as θ (1) φ(B) θ (1) Zt = β + at θ (B) φ(1) φ(1)

(A10)

so that a comparison of Equations (A9) and (A10) reveals that St =

φ(B) θ (1) Zt θ (B) φ(1)

(A11)

which allows St to be recursively estimated as St = −θ1 St−1 − · · · − θq St−q + (θ (1)/φ(1))(Zt + φ1 Zt−1 + · · · + φp Zt−p ), t > max(p, q)

(A12)

St = Zt ,

(A13)

t ≤ max(p, q)

A.2. Low-pass filtering techniques The HP filter is the solution to the problem of minimising

Nt2 + ξ

((St+1 − St ) − (St − St−1 ))2

(A14)

with respect to St and where ξ is a Lagrangean multiplier that has the interpretation of a smoothness parameter. The higher the value of ξ , the smoother the trend, so that in the limit, as ξ → ∞, St becomes a linear trend. The first-order conditions are 0 = −2(Zt − St ) + 2ξ((St − St−1 ) − (St−1 − St−2 )) − 4ξ((St+1 − St ) − (St − St−1 )) + 2ξ((St+2 − St+1 ) − (St+1 − St ))

(A15)

which may be written as Zt = ξ St−2 − 4ξ St−1 + (1 + 6ξ )St − 4ξ St+1 + ξ St+2

(A16)

so that Nt = Zt − St = ξ(St−2 − 4St−1 + 6St − 4St+1 + St+2 )

(A17)

Note that this expression cannot be used when t = 1, 2, T − 1, T so that, at these end-points of the sample, it is modified to N1 = ξ(S1 − 2S2 + S3 )

(A18a)

N2 = ξ(−2S1 + 5S2 − 4S3 + S4 )

(A18b)

NT −1 = ξ(ST −3 − 4ST −2 + 5ST −1 − 2ST )

(A18c)

NT = ξ(ST −2 − 2ST −1 + ST )

(A18d)


Int. J. Climatol. 26: 867–884 (2006)

882

T. C. MILLS

Defining Z = (Z1 , Z2 , . . . , ZT )T , S  1 −2 1 0 0  −2 5 −4 1 0   1 −4 6 −4 1   0 1 −4 6 −4  .  .  . = .  ..   0 ···   0 ···   0 ··· 0 ···

= (S1 , S2 , . . . , ST )T , N = (N1 , N2 , . . . , NT )T = Z − S and  0 0 ··· 0 0 0 ··· 0   0 0 ··· 0   1 0 ··· 0  ..   .  ..  .   0 1 −4 6 −4 1 0   0 0 1 −4 6 −4 1   0 0 0 1 −4 5 −2  0 0 0 0 1 −2 1

(A19)

the ‘system’ of Equations (A17) and (A18) can be written as N = ξ S

(A20)

Z = (I + ξ )S

(A21)

or

so that the vector of trend components is given by S = (I + ξ )−1 Z

(A22)

Equation (A16) can then be written as Zt = (ξ B 2 − 4ξ B + (1 + 6ξ ) − 4ξ B + ξ B 2 )St = (1 + ξ(1 − B)2 (1 − B −1 )2 )St

(A23)

so that the trend is St =

1 1 + ξ(1 − B)2 (1 − B −1 )2

Zt

(A24)

i.e. for an infinite sample, the HP filter estimates the trend as SˆtHP =

∞

wj Zt+j = wHP (B)Zt

(A25)

j =−∞

where wHP (B) =

1

(A26)

1 + ξ(1 − B)2 (1 − B −1 )2

For a finite sample, Equation (A22) can be used, but Equation (A25) is useful for investigating the frequency domain properties of the HP filter. Substituting B = e−iλ and B −1 = eiλ into Equation (A25) yields the frequency response function wHP (e−iλ ) =

1 1 + ξ(1 − e

−iλ 2

) (1 − e )

iλ 2

=

1 1 + 4ξ(1 − cos(λ))2

(A27)

using the trigonometric function eiλ + e−iλ = 2 cos(λ). Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)


883

Harvey and Jaeger (1993) show that the HP trend model Equation (A24) is the Weiner–Kolmogorov (WK) filter estimate of St in the ‘smooth trend plus white noise’ structural model, with the smoothing parameter being defined as the signal-to-noise ratio, i.e. (1 − B)2 St = νt

Nt = ut

ξ = σu2 /σν2

(A28)

Consider now an alternative unobserved component model of the form (1 − B)d St = (1 + B)f νt

Nt = (1 − B)f −d ut

(A29)

This has the WK trend filter estimate wB (B) =

1+ξ

1 (1 − B)(1 − B −1 )

f

(A30)

(1 + B)(1 + B −1 )

The filter wB (B) in Equation (A30) is known as the Butterworth (square-wave) filter. Although wB (B) is a two-sided infinite filter, Pollock (2000) presents an algorithm that enables the trend filter to be calculated in finite samples. The frequency response function of the Butterworth filter can be shown to be (see, e.g. Mills, 2003, pages 99–100) wB (e−iλ ) =

1

(A31)

1 + ξ(tan(λ/2))2f

As with the HP frequency response function of Equation (A27), higher values of ξ are associated with more abrupt cut-offs at lower frequencies, but in contrast to Equation (A27), the cut-offs using Equation (A31) are sharper. The smoothing parameter ξ in the HP and Butterworth filters may be chosen to correspond to a frequency cut-off of λc . Since the gain should equal 0.5 at λc , solving Equations (A27) and (A31) for this value gives ξHP = 1/4(1 − cos(λc ))2

(A32)

ξB = 1/(tan(λc /2))2f

(A33)

and

Given a choice of λc , and f in the case of the Butterworth filter, these expressions can be used to set values of the smoothing parameters in the two filters. REFERENCES Bacon DW, Watts DG. 1971. Estimating the transition between two intersecting straight lines. Biometrika 58: 525–534. Bell WR, Martin DEK. 2004. Computation of asymmetric signal extraction filters and mean squared error for ARIMA component models. Journal of Time Series Analysis 25: 603–625. Benestad RE. 2003. What can present climate models tell us about climate change? Climatic Change 59: 311–331. Beveridge S, Nelson CR. 1981. A new approach to decomposition of economic time series into permanent and transitory components with particular attention to measurement of the ‘business cycle’. Journal of Monetary Economics 7: 151–174. Brázdil R, Pfister C, Wanner H, von Storch H, Luterbacher J. 2005. Historical climatology in Europe – The state of the art. Climatic Change 70: 363–430. Craddock JM. 1957. An analysis of the slower temperature variations at Kew observatory by mean of exclusive band pass filters. Journal of The Royal Statistical Society. Series A 120: 387–397. Dickey DA, Fuller WA. 1979. Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association 74: 427–431. Essex C, McKitrick R. 2002. Taken by Storm. The Troubled Science, Policy and Politics of Global Warming. Key Porter Books: Toronto. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 867–884 (2006)

884

T. C. MILLS

Findlay DF, Monsell BC, Bell WR, Otto MC, Chen B-C. 1998. New capabilities and methods of the X-12-ARIMA seasonal-adjustment program. Journal of Business and Economic Statistics 16: 127–177. Gray AG, Thomson P. 2002. On a family of moving-average trend filters for the ends of a series. Journal of Forecasting 21: 125–149. Hall P, Tajvidi N. 2000. Nonparametric analysis of temporal trend when fitting parametric models to extreme-value data. Statistical Science 15: 153–167. Härdle W. 1990. Applied Nonparametric Regression. Cambridge University Press: Cambridge. Harvey AC. 1989. Forecasting, Structural Time Series Models and the Kalman filter. Cambridge University Press: Cambridge. Harvey AC. 2002. Trend analysis. In Encyclopedia of Environmetrics, El-Shaarawi AH, Piegorsch WW. (eds) Wiley: Chichester, 2243–2257. Harvey AC, Jaeger A. 1993. Detrending, stylized facts, and the business cycle. Journal of Applied Econometrics 8: 231–247. Harvey DI, Mills TC. 2001. Modelling global temperatures using cointegration and smooth transitions. Statistical Modelling 1: 143–159. Harvey DI, Mills TC. 2002. Unit roots and double smooth transitions. Journal of Applied Statistics 29: 675–683. Harvey DI, Mills TC. 2003. Modelling trends in Central England temperatures. Journal of Forecasting 22: 35–47. Hodrick RJ, Prescott EJ. 1997. Postwar US business cycles: an empirical investigation. Journal of Money, Credit and Banking 29: 1–16. Houghton JT, Ding Y, Griggs DJ, Nogeur M, van der Linden PJ, Dai X, Maskell K, Johnson CA. 2001. Intergovernmental Panel on Climate Change, Third Assessment Report: Climate Change 2001, The Scientific Basis, Cambridge University Press: Cambridge. Jones PD, Lister D. 2004. The development of monthly temperature series for Scotland and Northern Ireland. International Journal of Climatology 24: 569–590. Jones PD, Moberg A. 2003. Hemisphere and large-scale surface air temperature variation: An extensive revision and update to 2001. Journal of Climate 16: 206–223. Jones PD, New M, Parker DE, Martin S, Rigor IG. 1999. Surface air temperature and its changes over the past 150 years. Review of Geophysics 37: 173–199. Kendall MG. 1973. Time Series. Charles Griffin: London. Kendall MG, Stuart A, Ord JK. 1983. The Advanced Theory of Statistics. Vol. 3. Charles Griffin: London. Koopman SJ. 1993. Disturbance smoother for state space models. Biometrika 80: 117–126. Koopman SJ, Harvey AC, Doornik JA, Shephard N. 2000. STAMP: Structural Time Series Analyser, Modeller and Predictor. Timberlake Consultants: London. Leybourne SJ, Newbold P, Vougas D. 1998. Unit roots and smooth transitions. Journal of Time Series Analysis 19: 83–97. Mann ME. 2004. On smoothing potentially non-stationary climate time series. Geophysical Research Letters 31: L07214, DOI:10.1029/2004GL019569. Mann ME, Jones PD. 2003. Global surface temperatures over the past two millennia. Geophysical Research Letters 30: 1820, DOI:10.1029/2003GL017814. Mills TC. 2003. Modelling Trends and Cycles in Economic Time Series. Palgrave Macmillan: Basingstoke. Mills TC. 2004. Is the North Atlantic Oscillation a random walk? A comment with further results. International Journal of Climatology 24: 377–383. Mills TC. 2005. Modelling trends in England & Wales precipitation. Meteorological Applications 12: 169–176. Moberg A, Jones PD. 2005. Trends in indices for extremes in daily temperature and precipitation in central and Western Europe. International Journal of Climatology 25: 1149–1171. Newbold P, Vougas D. 1996. Beveridge-Nelson-type trends for I(2) and some seasonal models. Journal of Time Series Analysis 17: 151–169. Piˇsoft P, Kalvová J, Brázdil R. 2004. Cycles and trends in the Czech temperature series using wavelet transforms. International Journal of Climatology 24: 1661–1670. Pollock DSG. 1999. A Handbook of Time–Series Analysis, Signal Processing and Dynamics. Academic Press: San Diego, CA. Pollock DSG. 2000. Trend estimation and de-trending via rational square wave filters. Journal of Econometrics 99: 317–334. Quenneville B, Ladiray D, Lefran¸cois B. 2003. A note on Musgrave asymmetrical trend-cycle filters. International Journal of Forecasting 19: 727–734. Ruppert D, Wand MP, Carroll RJ. 2003. Semiparametric Regression. Cambridge University Press: Cambridge. Schmidlii J, Frei C. 2005. Trends of heavy precipitation and wet and dry spells in Switzerland during the 20th century. International Journal of Climatology 25: 753–771. Simonoff JS. 1996. Smoothing Methods in Statistics. Springer-Verlag: New York. Sollis R, Leybourne SJ, Newbold P. 1999. Unit roots and asymmetric smooth transitions. Journal of Time Series Analysis 20: 671–677. Soon WWH, Legates DR, Baliunas SL. 2004. Estimation and representation of long-term (>40 year) trends of Northern-Hemisphere gridded surface temperature. A note of caution. Geophysical Research Letters 31: L03209, DOI: 10.1029/2003GL019141. Stern DI, Kaufmann RK. 2000. Detecting a global warming signal in hemispheric temperature series: a structural time series analysis. Climatic Change 47: 411–438. Teräsvirta T. 1994. Specification, estimation, and evaluation of smooth transition autoregressive models. Journal of the American Statistical Association 89: 208–218. Visser H. 2004. Estimation and detection of flexible trends. Atmospheric Environment 38: 4135–4145. Whittaker ET. 1923. On a new method of graduations. Proceedings of the Edinburgh Mathematical Society 41: 63–75.


Int. J. Climatol. 26: 867–884 (2006)