Tendência STL de séries temporais usando R

27

Eu sou novo no R e na análise de séries temporais. Estou tentando encontrar a tendência de uma longa série temporal de temperatura diária (40 anos) e tentei aproximações diferentes. O primeiro é apenas uma regressão linear simples e o segundo é a Decomposição Sazonal de Séries Temporais por Loess.

Neste último, parece que o componente sazonal é maior que a tendência. Mas, como quantifico a tendência? Gostaria apenas de um número dizendo o quão forte é essa tendência.

     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

insira a descrição da imagem aqui

pacomet
fonte

Respostas:

20

Eu não me incomodaria com stl()isso - a largura de banda do lowess mais suave usada para extrair a tendência é muito, muito pequena, resultando nas flutuações de pequena escala que você vê. Eu usaria um modelo aditivo. Aqui está um exemplo usando o código de dados e modelo do livro de Simon Wood sobre 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")

dados de temperatura cairo

Ajuste um modelo com tendência e componentes sazonais --- aviso de que isso é lento:

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

O modelo ajustado é assim:

> 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

e podemos visualizar a tendência e os termos sazonais via

plot(mod$gam, pages = 1)

Cairo tendência ajustada e sazonal

e se quisermos traçar a tendência nos dados observados, podemos fazer isso com previsão 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)

Tendência do Cairo

Ou o mesmo para o modelo real:

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)

Modelo montado no Cairo

Este é apenas um exemplo, e uma análise mais aprofundada pode ter que lidar com o fato de que existem alguns dados ausentes, mas o acima deve ser um bom ponto de partida.

Quanto ao seu argumento sobre como quantificar a tendência - bem, isso é um problema, porque a tendência não é linear, nem na sua stl()versão nem na versão GAM que mostro. Se fosse, você poderia indicar a taxa de variação (inclinação). Se você deseja saber até que ponto a tendência estimada mudou ao longo do período de amostragem, podemos usar os dados contidos prede calcular a diferença entre o início e o final da série apenas no componente de tendência :

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

então as temperaturas são, em média, 1,76 graus mais quentes do que no início do registro.

Restabelecer Monica - G. Simpson
fonte
Olhando para o gráfico, acho que pode haver alguma confusão entre Fahrenheit e Celsius.
Henry
Bem visto - estou fazendo algo semelhante há alguns meses e os dados estão no grau C. Era força do hábito!
Reintegrar Monica - G. Simpson
Obrigado Gavin, uma resposta muito agradável e compreensível. Vou tentar suas sugestões. É uma boa idéia plotar o componente stl () trend e fazer uma regressão linear?
pacomet
1
@ pacomet - não, não realmente, a menos que você ajuste um modelo que seja responsável pela autocorrelação nos resíduos, como fiz acima. Você pode usar o GLS para isso ( gls()no pacote nlme). Mas, como mostra acima, no Cairo, e o STL sugere para seus dados, a tendência não é linear. Como tal, uma tendência linear não seria apropriada - pois não descreve os dados corretamente. Você precisa experimentá-lo com seus dados, mas um AM como eu mostro seria degradado para uma tendência linear se isso se ajustasse melhor aos dados.
Reintegrar Monica - G. Simpson
1
@ andreas-h eu não faria isso; a tendência STL está super ajustada. Ajuste o GAM com a estrutura AR () e interprete a tendência. Isso fornecerá um modelo de regressão adequado que será muito mais útil para você.
Reponha Monica - G. Simpson
4

Gavin forneceu uma resposta muito completa, mas para uma solução mais simples e rápida, recomendo definir o parâmetro da função stl t.window para um valor que seja múltiplo da frequência dos dados ts . Eu usaria a periodicidade inferida de interesse (por exemplo, um valor de 3660 para tendências decadais com dados de resolução diurnos). Você também pode estar interessado no pacote stl2 descrito na dissertação do autor . Eu apliquei o método de Gavin aos meus próprios dados e também é muito eficaz.

Adam Erickson
fonte