Tendencia STL de series de tiempo usando R

Soy nuevo en R y en el análisis de series de tiempo. Estoy tratando de encontrar la tendencia de una serie temporal de temperatura diaria larga (40 años) e intenté diferentes aproximaciones. La primera es solo una regresión lineal simple y la segunda es la descomposición estacional de series temporales de Loess.

En este último parece que el componente estacional es mayor que la tendencia. Pero, ¿cómo cuantifico la tendencia? Me gustaría solo un número que diga qué tan fuerte es esa tendencia.

     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

ingrese la descripción de la imagen aquí

r time-series trend pacomet
fuente

No me molestaría con stl()esto: el ancho de banda para el lowess más suave utilizado para extraer la tendencia es muy, muy pequeño, lo que resulta en las fluctuaciones de pequeña escala que ves. Usaría un modelo aditivo. Aquí hay un ejemplo que usa datos y código de modelo del libro de Simon Wood sobre 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")

datos de temperatura del cairo

Ajuste un modelo con tendencia y componentes estacionales --- advertencia esto es 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)))

El modelo ajustado se ve así:

> 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

y podemos visualizar la tendencia y los términos estacionales a través de

plot(mod$gam, pages = 1)

El Cairo equipado tendencia y estacional

y si queremos trazar la tendencia en los datos observados podemos hacerlo con predicción a través de:

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)

Tendencia ajustada de El Cairo

O lo mismo para el 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 ajustado Cairo

Este es solo un ejemplo, y un análisis más profundo podría tener que lidiar con el hecho de que faltan algunos datos, pero lo anterior debería ser un buen punto de partida.

En cuanto a su punto sobre cómo cuantificar la tendencia, bueno, eso es un problema, porque la tendencia no es lineal, ni en su stl()versión ni en la versión GAM que muestro. Si lo fuera, podría dar la tasa de cambio (pendiente). Si desea saber cuánto ha cambiado la tendencia estimada durante el período de muestreo, entonces podemos usar los datos contenidos predy calcular la diferencia entre el inicio y el final de la serie solo en el componente de tendencia :

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

entonces las temperaturas son, en promedio, 1.76 grados más cálidas que al comienzo del registro.

Restablece a Mónica - G. Simpson
fuente

Al mirar el gráfico, creo que puede haber cierta confusión entre Fahrenheit y Celsius.

Henry

Bien visto: he estado haciendo algo similar durante unos meses y los datos están en grado C. ¡Era la fuerza del hábito!

Restablece a Monica - G. Simpson el

Gracias Gavin, una respuesta muy agradable y comprensible. Probaré tus sugerencias. ¿Es una buena idea trazar el componente de tendencia stl () y hacer una regresión lineal?

pacomet

@pacomet: no, en realidad no, a menos que se ajuste a un modelo que tenga en cuenta la autocorrelación en los residuos como lo hice anteriormente. Puede usar GLS para eso ( gls()en el paquete nlme). Pero como lo anterior muestra para El Cairo, y el STL sugiere para sus datos, la tendencia no es lineal. Como tal, una tendencia lineal no sería apropiada, ya que no describe los datos correctamente. Debe probarlo con sus datos, pero una AM como la muestro se degradaría a una tendencia lineal si se ajustara mejor a los datos.

Reinstale a Monica - G. Simpson el

@ andreas-h no haría eso; La tendencia STL está demasiado ajustada. Ajuste el GAM con la estructura AR () e interprete la tendencia. Eso le dará un modelo de regresión adecuado que será mucho más útil para usted.

Restablecer a Monica - G. Simpson

Tendencia STL de series de tiempo usando R

Respuestas: