STL тренд временных рядов с использованием R

27

Я новичок в R и для анализа временных рядов. Я пытаюсь найти тренд длинного (40 лет) дневного временного ряда температуры и пробовал разные приближения. Первый - это простая линейная регрессия, а второй - сезонная декомпозиция временных рядов по Лесс.

В последнем случае оказывается, что сезонная составляющая больше, чем тренд. Но как мне определить тренд? Я хотел бы просто указать, насколько сильна эта тенденция.

     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

введите описание изображения здесь

pacomet
источник

Ответы:

20

Я бы не стал беспокоиться stl()об этом - пропускная способность низкочастотного сглаживателя, используемого для извлечения тренда, очень мала, что приводит к небольшим колебаниям, которые вы видите. Я бы использовал аддитивную модель. Вот пример использования данных и кода модели из книги Саймона Вуда о 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")

Каир данные о температуре

Установите модель с трендовыми и сезонными компонентами - предупреждая, что это медленно:

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

Подогнанная модель выглядит так:

> 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

и мы можем визуализировать тенденцию и сезонные условия с помощью

plot(mod$gam, pages = 1)

Каир соответствует тренду и сезонности

и если мы хотим построить тренд на наблюдаемых данных, мы можем сделать это с помощью прогноза:

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)

Каир соответствует тренду

Или то же самое для фактической модели:

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)

Каир подогнанная модель

Это всего лишь пример, и более глубокому анализу, возможно, придется иметь дело с тем фактом, что существует несколько пропущенных данных, но вышеизложенное должно стать хорошей отправной точкой.

Что касается вашей точки зрения о том, как количественно оценить тренд - это проблема, потому что тренд не является линейным, ни в вашей stl()версии, ни в версии GAM, которую я показываю. Если бы это было так, вы могли бы указать скорость изменения (уклон). Если вы хотите узнать, насколько изменился предполагаемый тренд за период выборки, мы можем использовать данные, содержащиеся в нем, predи вычислить разницу между началом и концом ряда только в компоненте тренда :

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

таким образом, температура в среднем на 1,76 градуса выше, чем в начале записи.

Восстановить Монику - Дж. Симпсон
источник
Глядя на график, я думаю, что может быть некоторая путаница между Фаренгейтом и Цельсией.
Генри
Хорошо заметили - я делал что-то подобное в течение нескольких месяцев, и данные в степени C. Была сила привычки!
Восстановить Монику - Дж. Симпсон
Спасибо Гэвин, очень хороший и понятный ответ. Я попробую ваши предложения. Хорошая идея построить график трендового компонента stl () и сделать линейную регрессию?
pacomet
1
@pacomet - нет, не совсем, если только вы не подходите к модели, которая учитывает автокорреляцию в остатках, как я делал выше. Вы можете использовать GLS для этого ( gls()в пакете NLME). Но, как показано выше для Каира, и STL предлагает для ваших данных, тенденция не является линейной. Таким образом, линейный тренд не подходит, поскольку он не описывает данные должным образом. Вы должны попробовать это на своих данных, но AM, как я показываю, будет ухудшаться до линейного тренда, если он лучше всего соответствует данным.
Восстановить Монику - Дж. Симпсон
1
@ andreas-h я бы этого не делал; тенденция STL переоценена. Установите GAM со структурой AR () и интерпретируйте тренд. Это даст правильную модель регрессии, которая будет гораздо более полезной для вас.
Восстановить Монику - Г. Симпсон
4

Гэвин дал очень подробный ответ, но для более простого и быстрого решения я рекомендую установить для параметра stl function t.window значение, кратное частоте данных ts . Я бы использовал предполагаемую периодичность интереса (например, значение 3660 для десятилетних трендов с данными суточного разрешения). Вас также может заинтересовать пакет stl2 , описанный в диссертации автора . Я применил метод Гэвина к своим собственным данным, и он тоже очень эффективен.

Адам Эриксон
источник