Trend STL szeregów czasowych przy użyciu R.

Jestem nowy w R i analizie szeregów czasowych. Próbuję znaleźć trend w długim (40 lat) dziennym szeregu czasowym temperatur i próbowałem różnych przybliżeń. Pierwszy to po prostu prosta regresja liniowa, a drugi to Sezonowy rozkład szeregów czasowych według Loessa.

W tym ostatnim wydaje się, że składnik sezonowy jest większy niż trend. Ale jak mam zmierzyć trend? Chciałbym tylko powiedzieć, jak silny jest ten 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

wprowadź opis zdjęcia tutaj

r time-series trend pacomet
źródło

Nie zawracałbym stl()sobie tym głowy - szerokość pasma dla bezużytecznej wygładzarki wykorzystywanej do wydobywania trendu jest daleko, daleko, aż do małych, co powoduje fluktuacje na małą skalę, które widzisz. Użyłbym modelu addytywnego. Oto przykład wykorzystujący dane i kod modelu z książki Simona Wooda na temat 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")

dane temperatury Cairo

Dopasuj model do trendu i komponentów sezonowych --- ostrzeżenie, że jest to powolne:

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

Dopasowany model wygląda następująco:

> 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

i możemy wizualizować trendy i warunki sezonowe za pośrednictwem

plot(mod$gam, pages = 1)

Kair dopasował trend i sezonowość

a jeśli chcemy wykreślić trend na obserwowanych danych, możemy to zrobić z prognozą poprzez:

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)

Trend dopasowany do Kairu

Lub to samo dla rzeczywistego modelu:

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)

Model dopasowany do Kairu

To tylko przykład, a dogłębniejsza analiza może wymagać uwzględnienia kilku brakujących danych, ale powyższe powinno być dobrym punktem wyjścia.

Jeśli chodzi o twoje zdanie na temat kwantyfikacji trendu - to jest problem, ponieważ trend nie jest liniowy, ani w twojej stl()wersji, ani w wersji GAM, którą pokazuję. Jeśli tak, możesz podać tempo zmian (nachylenie). Jeśli chcesz wiedzieć, o ile zmienił się szacowany trend w okresie próbkowania, możemy wykorzystać zawarte w nim dane predi obliczyć różnicę między początkiem i końcem serii tylko w komponencie trendu :

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

więc temperatury są średnio o 1,76 stopnia wyższe niż na początku zapisu.

Przywróć Monikę - G. Simpson
źródło

Patrząc na mapę, myślę, że może być pewne zamieszanie między Fahrenheitem a Celsjuszem.

Henry

Dobrze zauważony - robiłem coś podobnego od kilku miesięcy, a dane są w stopniu C. Była siłą nawyku!

Przywróć Monikę - G. Simpson

Dzięki Gavin, bardzo miła i zrozumiała odpowiedź. Spróbuję twoich sugestii. Czy dobrym pomysłem jest wykreślenie komponentu trendu stl () i wykonanie regresji liniowej?

pacomet

@pacomet - nie, nie do końca, chyba że pasujesz do modelu uwzględniającego autokorelację w resztkach, tak jak ja powyżej. Możesz do tego użyć GLS ( gls()w pakiecie nlme). Ale jak pokazano powyżej dla Kairu i STL sugeruje dla twoich danych, trend nie jest liniowy. Jako taki trend liniowy nie byłby odpowiedni, ponieważ nie opisuje poprawnie danych. Musisz wypróbować to na swoich danych, ale AM, taki jak ja, pokazałby się jako trend liniowy, gdyby najlepiej pasował do danych.

Przywróć Monikę - G. Simpson

@ andreas-h Nie zrobiłbym tego; tendencja STL jest już przesadzona. Dopasuj GAM do struktury AR () i zinterpretuj trend. To da odpowiedni model regresji, który będzie ci znacznie bardziej przydatny.

Przywróć Monikę - G. Simpson

Trend STL szeregów czasowych przy użyciu R.

Odpowiedzi: