STL tendencia de las series de tiempo mediante R

Question

Soy nuevo en R y análisis de series de tiempo. Estoy tratando de encontrar la tendencia de largo (40 años) la temperatura diaria de series de tiempo y trató de aproximaciones diferentes. Primero es sólo una de regresión lineal simple y la segunda es la Temporada de Descomposición de Series de Tiempo por el Loess.

En este último parece que el componente estacional es mayor que el de la tendencia. Pero, ¿cómo cuantificar la tendencia? Me gustaría simplemente un número que le dice 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

Answer 1

Que no me molestaría en stl() - el ancho de banda para la lowess smoother utiliza para extraer la tendencia es lejos, lejos, a los pequeños resultante de las fluctuaciones en pequeña escala que ver. Me gustaría utilizar un modelo aditivo. Aquí hay un ejemplo usando los datos y el código de modelo de Simon de Madera, el libro de 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")

Ajuste de un modelo con tendencia y estacionalidad --- advertencia este es lenta:

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 parece a esto:

> 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 estacional de los términos a través de

plot(mod$gam, pages = 1)

y si queremos representar la tendencia en los datos observados podemos hacerlo con la 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)

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

Este es sólo un ejemplo, y un análisis más profundo podría tener que lidiar con el hecho de que hay un par de datos que faltan, pero la de arriba debe ser un buen punto de partida.

En cuanto a su punto acerca de cómo cuantificar la tendencia - bueno, eso es un problema, porque la tendencia no es lineal, ni en su stl() versión ni la GAM versión muestro. Si lo fuera, usted podría dar a la tasa de cambio (pendiente). Si quieres saber por cuánto se ha estimado observa un cambio de tendencia durante el período de muestreo, entonces podemos usar los datos contenidos en pred y calcular la diferencia entre el inicio y el final de la serie en la tendencia sólo componentes:

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

por lo que las temperaturas son, en promedio, 1.76 grados más caliente que en el inicio de la grabación.

Answer 2

Gavin se proporciona una respuesta detallada, pero de una forma más sencilla y más rápida solución, le recomendamos que ajuste la stl función de t.ventana de parámetro a un valor que es un múltiplo de la frecuencia de la ts de datos. Me gustaría utilizar el inferirse periodicidad de interés (por ejemplo, un valor de 3660 para decadal de las tendencias diurna resolución de datos). Usted también podría estar interesado en el stl2 paquete descrito en el autor de la tesis doctoral. He aplicado Gavin del método para mis propios datos y es muy eficaz también.