Tendance STL des séries chronologiques utilisant R

27

Je suis novice en R et en analyse de séries chronologiques. J'essaie de trouver la tendance d'une longue série de températures quotidiennes (40 ans) et j'ai essayé différentes approximations. Le premier n'est qu'une simple régression linéaire et le second est la décomposition saisonnière des séries chronologiques par Loess.

Dans ce dernier, il apparaît que la composante saisonnière est supérieure à la tendance. Mais comment quantifier la tendance? Je voudrais juste un chiffre indiquant à quel point cette tendance est forte.

     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

entrez la description de l'image ici

pacomet
la source

Réponses:

20

Je ne m'embêterais pas stl()pour cela - la bande passante pour le lisseur lowess utilisé pour extraire la tendance est loin, très petite, ce qui entraîne les petites fluctuations d'échelle que vous voyez. J'utiliserais un modèle additif. Voici un exemple utilisant les données et le code modèle du livre de Simon Wood sur les GAM:

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")

données de température du Caire

Ajustez un modèle avec des composants tendance et saisonniers --- attention, c'est lent:

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)))

Le modèle ajusté ressemble à ceci:

> 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

et nous pouvons visualiser la tendance et les conditions saisonnières via

plot(mod$gam, pages = 1)

Caire tendance ajustée et saisonnière

et si nous voulons tracer la tendance sur les données observées, nous pouvons le faire avec la prédiction 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)

Tendance ajustée du Caire

Ou la même chose pour le modèle réel:

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)

Modèle ajusté Cairo

Ce n'est qu'un exemple, et une analyse plus approfondie pourrait devoir traiter du fait qu'il manque quelques données, mais ce qui précède devrait être un bon point de départ.

Quant à votre point sur la façon de quantifier la tendance - eh bien c'est un problème, car la tendance n'est pas linéaire, ni dans votre stl()version ni dans la version GAM que je montre. Si c'était le cas, vous pourriez donner le taux de changement (pente). Si vous voulez savoir dans quelle mesure la tendance estimée a changé au cours de la période d'échantillonnage, nous pouvons utiliser les données contenues dans predet calculer la différence entre le début et la fin de la série dans la composante tendance uniquement:

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

les températures sont donc, en moyenne, 1,76 degrés plus chaudes qu'au début du record.

Réintégrer Monica - G. Simpson
la source
En regardant le tableau, je pense qu'il peut y avoir une certaine confusion entre Fahrenheit et Celsius.
Henry
Bien repéré - je fais quelque chose de similaire depuis quelques mois et les données sont en degré C. C'était une force d'habitude!
Reinstate Monica - G. Simpson
Merci Gavin, une réponse très agréable et compréhensible. J'essaierai vos suggestions. Est-ce une bonne idée de tracer la composante de tendance stl () et de faire une régression linéaire?
pacomet
1
@pacomet - non, pas vraiment, sauf si vous ajustez un modèle qui tient compte de l'autocorrélation dans les résidus comme je l'ai fait ci-dessus. Vous pouvez utiliser GLS pour cela ( gls()dans le package nlme). Mais comme le montre ce qui précède pour le Caire, et la STL le suggère pour vos données, la tendance n'est pas linéaire. En tant que tel, une tendance linéaire ne serait pas appropriée - car elle ne décrit pas correctement les données. Vous devez l'essayer sur vos données, mais un AM comme je le montre se dégraderait en une tendance linéaire si cela correspondait le mieux aux données.
Reinstate Monica - G. Simpson
1
@ andreas-h je ne ferais pas ça; la tendance STL est trop ajustée. Ajustez le GAM avec la structure AR () et interprétez la tendance. Cela vous donnera un modèle de régression approprié qui vous sera beaucoup plus utile.
Rétablir Monica - G. Simpson
4

Gavin a fourni une réponse très approfondie, mais pour une solution plus simple et plus rapide, je recommande de définir le paramètre stl function t.window sur une valeur qui est un multiple de la fréquence des données ts . J'utiliserais la périodicité présumée d'intérêt (par exemple, une valeur de 3660 pour les tendances décennales avec des données de résolution diurne). Vous pouvez également être intéressé par le package stl2 décrit dans la dissertation de l'auteur . J'ai appliqué la méthode de Gavin à mes propres données et elle est également très efficace.

Adam Erickson
la source