Tendência STL de séries temporais usando R

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

r time-series trend pacomet
fonte

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

@ 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

@ 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

Tendência STL de séries temporais usando R

Respostas: