Tren time series STL menggunakan R

Saya baru mengenal R dan analisis deret waktu. Saya mencoba menemukan tren rangkaian waktu suhu harian yang panjang (40 tahun) dan mencoba perkiraan yang berbeda. Yang pertama hanyalah regresi linier sederhana dan yang kedua adalah Dekomposisi Musiman dari Time Series oleh Loess.

Dalam yang terakhir tampak bahwa komponen musiman lebih besar daripada tren. Tapi, bagaimana saya mengukur tren? Saya hanya ingin angka yang mengatakan seberapa kuat tren itu.

     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

masukkan deskripsi gambar di sini

r time-series trend pacomet
sumber

Saya tidak akan repot dengan stl()ini - bandwidth untuk lowess yang lebih halus yang digunakan untuk mengekstraksi tren jauh, jauh, menjadi kecil yang menghasilkan fluktuasi skala kecil yang Anda lihat. Saya akan menggunakan model tambahan. Berikut ini adalah contoh menggunakan data dan kode model dari buku Simon Wood tentang 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")

data suhu cairo

Pas dengan model dengan komponen tren dan musiman --- peringatan ini lambat:

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

Model yang pas terlihat seperti ini:

> 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

dan kami dapat memvisualisasikan tren dan ketentuan musiman melalui

plot(mod$gam, pages = 1)

Tren Kairo pas dan musiman

dan jika kita ingin merencanakan tren pada data yang diamati, kita dapat melakukannya dengan prediksi melalui:

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)

Kairo cocok dengan tren

Atau sama untuk model aktual:

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 Kairo pas

Ini hanya contoh, dan analisis yang lebih mendalam mungkin harus berurusan dengan fakta bahwa ada beberapa data yang hilang, tetapi di atas harus menjadi titik awal yang baik.

Mengenai poin Anda tentang bagaimana mengukur tren - yah itu masalah, karena tren itu tidak linier, baik dalam stl()versi Anda maupun versi GAM yang saya tunjukkan. Jika ya, Anda bisa memberikan tingkat perubahan (kemiringan). Jika Anda ingin tahu seberapa besar kecenderungan tren berubah selama periode pengambilan sampel, maka kita dapat menggunakan data yang terkandung di dalamnya preddan menghitung perbedaan antara awal dan akhir seri hanya dalam komponen tren :

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

jadi suhu rata-rata 1,76 derajat lebih hangat daripada di awal catatan.

Pasang kembali Monica - G. Simpson
sumber

Melihat grafik, saya pikir mungkin ada beberapa kebingungan antara Fahrenheit dan Celsius.

Henry

Terlihat dengan baik - Saya telah melakukan hal yang serupa selama beberapa bulan dan datanya dalam derajat C. Sudah menjadi kebiasaan!

Pasang kembali Monica - G. Simpson

Terima kasih Gavin, jawaban yang sangat bagus dan mudah dimengerti. Saya akan mencoba saran Anda. Apakah ide yang baik untuk merencanakan komponen tren stl () dan membuat regresi linier?

pacomet

@ometomet - tidak, tidak terlalu, kecuali jika Anda cocok dengan model yang memperhitungkan autokorelasi dalam residu seperti yang saya lakukan di atas. Anda dapat menggunakan GLS untuk itu ( gls()dalam paket nlme). Tetapi seperti yang ditunjukkan di atas untuk Kairo, dan STL menyarankan untuk data Anda, trennya tidak linear. Dengan demikian, tren linier tidak akan sesuai - karena tidak menggambarkan data dengan benar. Anda perlu mencobanya pada data Anda, tetapi AM seperti yang saya tunjukkan akan menurunkan tren linier jika itu cocok dengan data terbaik.

Pasang kembali Monica - G. Simpson

@ andreas-h Saya tidak akan melakukan itu; tren STL sudah terlalu pas. Pasangkan GAM dengan struktur AR () dan interpretasikan tren. Itu akan memberikan model regresi yang tepat yang akan jauh lebih bermanfaat bagi Anda.

Pasang kembali Monica - G. Simpson

Tren time series STL menggunakan R

Jawaban: