27
$\begingroup$

I am new to R and to time series analysis. I am trying to find the trend of a long (40 years) daily temperature time series and tried to different approximations. First one is just a simple linear regression and second one is Seasonal Decomposition of Time Series by Loess.

In the latter it appears that the seasonal component is greater than the trend. But, how do I quantify the trend? I would like just a number telling how strong is that trend.

     Call:  stl(x = tsdata, s.window = "periodic")
     Time.series components:
        seasonal                trend            remainder               
Min.   :-8.482470191   Min.   :20.76670   Min.   :-11.863290365      
1st Qu.:-5.799037090   1st Qu.:22.17939   1st Qu.: -1.661246674 
Median :-0.756729578   Median :22.56694   Median :  0.026579468      
Mean   :-0.005442784   Mean   :22.53063   Mean   : -0.003716813 
3rd Qu.:5.695720249    3rd Qu.:22.91756   3rd Qu.:  1.700826647    
Max.   :9.919315613    Max.   :24.98834   Max.   : 12.305103891   

 IQR:
         STL.seasonal STL.trend STL.remainder data   
         11.4948       0.7382    3.3621       10.8051
       % 106.4          6.8      31.1         100.0  
     Weights: all == 1
     Other components: List of 5   
$ win  : Named num [1:3] 153411 549 365  
$ deg  : Named int [1:3] 0 1 1   
$ jump : Named num [1:3] 15342 55 37  
$ inner: int 2  
$ outer: int 0

enter image description here

$\endgroup$
20
$\begingroup$

I wouldn't bother with stl() for this - the bandwidth for the lowess smoother used to extract the trend is far, far, to small resulting in the small scale fluctuations you see. I would use an additive model. Here is an example using data and model code from Simon Wood's book on GAMs:

require(mgcv)
require(gamair)
data(cairo)
cairo2 <- within(cairo, Date <- as.Date(paste(year, month, day.of.month, 
                                              sep = "-")))
plot(temp ~ Date, data = cairo2, type = "l")

cairo temperature data

Fit a model with trend and seasonal components --- warning this is slow:

mod <- gamm(temp ~ s(day.of.year, bs = "cc") + s(time, bs = "cr"),
            data = cairo2, method = "REML",
            correlation = corAR1(form = ~ 1 | year),
            knots = list(day.of.year = c(0, 366)))

The fitted model looks like this:

> summary(mod$gam)

Family: gaussian 
Link function: identity 

Formula:
temp ~ s(day.of.year, bs = "cc") + s(time, bs = "cr")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  71.6603     0.1523   470.7   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Approximate significance of smooth terms:
                 edf Ref.df       F p-value    
s(day.of.year) 7.092  7.092 555.407 < 2e-16 ***
s(time)        1.383  1.383   7.035 0.00345 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

R-sq.(adj) =  0.848  Scale est. = 16.572    n = 3780

and we can visualise the trend and seasonal terms via

plot(mod$gam, pages = 1)

Cairo fitted trend and seasonal

and if we want to plot the trend on the observed data we can do that with prediction via:

pred <- predict(mod$gam, newdata = cairo2, type = "terms")
ptemp <- attr(pred, "constant") + pred[,2]
plot(temp ~ Date, data = cairo2, type = "l",
     xlab = "year",
     ylab = expression(Temperature ~ (degree*F)))
lines(ptemp ~ Date, data = cairo2, col = "red", lwd = 2)

Cairo fitted trend

Or the same for the actual model:

pred2 <- predict(mod$gam, newdata = cairo2)
plot(temp ~ Date, data = cairo2, type = "l",
     xlab = "year",
     ylab = expression(Temperature ~ (degree*F)))
lines(pred2 ~ Date, data = cairo2, col = "red", lwd = 2)

Cairo fitted model

This is just an example, and a more in-depth analysis might have to deal with the fact that there are a few missing data, but the above should be a good starting point.

As to your point about how to quantify the trend - well that is a problem, because the trend is not linear, neither in your stl() version nor the GAM version I show. If it were, you could give the rate of change (slope). If you want to know by how much has the estimated trend changed over the period of sampling, then we can use the data contained in pred and compute the difference between the start and the end of the series in the trend component only:

> tail(pred[,2], 1) - head(pred[,2], 1)
    3794 
1.756163

so temperatures are, on average, 1.76 degrees warmer than at the start of the record.

$\endgroup$
  • $\begingroup$ Looking at the chart, I think there may be some confusion between Fahrenheit and Celsius. $\endgroup$ – Henry Apr 13 '11 at 15:36
  • $\begingroup$ Well spotted - I've been doing something similar for a few months and the data are in degree C. Was force of habit! $\endgroup$ – Gavin Simpson Apr 13 '11 at 15:40
  • $\begingroup$ Thanks Gavin, a very nice and understandable answer. I'll try your suggestions. Is it a good idea to plot the stl() trend component and make a linear regression? $\endgroup$ – pacomet Apr 14 '11 at 7:59
  • 1
    $\begingroup$ @pacomet - no, not really, unless you fit a model that accounts for autocorrelation in the residuals as I did above. You could use GLS for that (gls() in package nlme). But as the above shows for Cairo, and the STL suggests for your data, the trend isn't linear. As such, a linear trend would not be appropriate - as it doesn't describe the data properly. You need to try it on your data, but an AM like I show would degrade to a linear trend if that fitted the data best. $\endgroup$ – Gavin Simpson Apr 14 '11 at 8:07
  • 1
    $\begingroup$ @andreas-h I wouldn't do that; the STL trend is over fitted. Fit the GAM with the AR() structure and interpret the trend. That will give a proper regression model which will be far more useful to you. $\endgroup$ – Gavin Simpson May 17 '12 at 13:56
4
$\begingroup$

Gavin provided a very thorough answer, but for a simpler and faster solution, I recommend setting the stl function t.window parameter to a value that is a multiple of the frequency of the ts data. I would use the inferred periodicity of interest (e.g., a value of 3660 for decadal trends with diurnal resolution data). You may also be interested in the stl2 package described in the author's dissertation. I have applied Gavin's method to my own data and it is very effective too.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.