1. Что такое периодограмма?

Пусть \(X_N = (x_1, \ldots, x_N)\) — ряд длины \(N\). \[ x_n = c_0 + \sum_{k = 1}^{[N/2]} \left( c_k\cos(2\pi nk/N) + s_k \sin(2\pi nk/N) \right), \] Разложение Фурье ряда \(X_N\). Здесь \(0 \le n < N\) и \(s_{N/2} = 0\) для четного \(N\).

Периодограммой ряда \(X_N\) называется функция \(П_x^N\), которая задается формулой \[\begin{equation*} П_x^N(k/N) = \frac{N}{2} \begin{cases} 2c_0^2 &\text{для k = 0},\\ c_k^2 + s_k^2 &\text{для 0<k<N/2},\\ 2c_{N/2}^2 &\text{для k = N/2}. \end{cases} \end{equation*}\]

Норма ряда имеет вид \[ || X_N ||^2 = \sum_{k = 0}^{[N/2]} П_x^N(k/N) \] То есть значение периодограммы в точке \(k/N\) описывает вклад гармонической компоненты с частотой \(\omega = k/N\). По периодограмме можно понять, какие частоты есть в разложении Фурье данного ряда.

2. Что такое тренд? Периодограмма тренда.

Тренд - основная тенденция изменения временного ряда. Если рассматривать ряд, как случайный процесс, то тренд можно определить как математическое ожидание процесса. Если ряд рассматривать как функцию (детерминированную), то тренд можно аппроксимировать некоторой функцией, например полиномом.

Периодограмма тренда для смоделированных данных.

x <- 1:200
y <- exp(0.1*x)

sp <- spec.pgram(y,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

3. Что такое периодическая компонента? Ее периодограмма.

Сезонной компонентой называют периодически повторяющуюся компоненту.

Пример периодограммы периодической компоненты для смоделированных данных.

x <- 1:200
y <- cos(2*pi*(1/12)*x) + cos(2*pi*(2/12)*x) + cos(2*pi*(3/12)*x)+ cos(2*pi*(4/12)*x) + cos(2*pi*(5/12)*x)
str <- c("0","1/12","2/12","3/12","4/12","5/12","6/12")
sp <- spec.pgram(y,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0, plot = FALSE)
{plot(sp,log = "no",xaxt = "n")
      axis(1,at = seq(0,0.5,by = 1/12),las = 2,labels=sprintf("%s",str))}

4. Что такое шум? Как выглядит его периодограмма?

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

Пример периодограммы белого шума для смоделированных данных.

y <- rnorm(1000)

spec.pgram(y,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

5. Периодограмма как оценка спектральной плотности. Распределение значений. Сглаживание периодограммы

Если процесс стационарный, то по периодограмме можно строить оценки спектральной плотности ряда. Так как на практике данные зашумлены и не очень важно точно определять частоты основных функций синусов и косинусов, имеет смысл сгладить периодограмму, для того чтобы выделить частоты с наибольшим вкладом в периодику ряда, при этом уменьшив значение случайных пиков. Существует много различных подходов к задаче сглаживании периодограммы. В R в функции spec.pgram реализовано сглаживание с окном Даниэля. Окно Даниэля означает простое (с равными весами) сглаживание скользящим средним значений периодограммы.

Пример сглаживания периодограммы на модельном ряде с шумом.

x <- 1:200
y <- cos(2*pi*(1/12)*x) + cos(2*pi*(2/12)*x) + cos(2*pi*(3/12)*x)+ cos(2*pi*(4/12)*x) + cos(2*pi*(5/12)*x) + rnorm(200)
sp <- spec.pgram(y,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0, plot = FALSE)
{plot(sp,log = "no",xaxt = "n", main = "Несглаженная периодограмма")
      axis(1,at = seq(0,0.5,by = 1/12),las = 2,labels=sprintf("%s",str))}

sp <- spec.pgram(y,spans = c(3,3), detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0, plot = FALSE)
{plot(sp,log = "no",xaxt = "n",  main = "Сглаженная периодограмма")
      axis(1,at = seq(0,0.5,by = 1/12),las = 2,labels=sprintf("%s",str))}

6. Чем отличается сглаживание от выделения тренда?

При сглаживании ряда мы получаем медленно меняющуюся компоненту. О тренде можно говорить как о случайном или неслучайном. Если под трендом понимаем детерменированную функцию, то эта функция будет трендом (или нет).

7. Что такое линейный фильтр, импульсная характеристика? Причинный фильтр, FIR.

Пусть \(x_n\) - последовательность.

Линейным фильтром называется функция \[(\Phi(X))_j = \sum_{i = -\infty}^\infty h_i x_{j-1},\] здесь \({h_i}\) - импульсная характеристика.

FIR - фильтр с конечной импульсной характеристикой.

Причинный фильтр определяется выражением \[(\Phi(X))_j = \sum_{i = 0}^r h_i x_{j-1}.\]

8. Характеристики фильтра через его воздействие на \(\cos(2\pi \omega n)\) (или на комплексную экспоненту) – АЧХ, ФЧХ.

\(H_\Phi(z) = \sum_i h_i z^{-i}\) - передаточная функция.

Пусть \(z= \exp (i2\pi \omega)\), тогда \(|H_\Phi(\exp (i2\pi \omega))| = A_\Phi(\omega)\) - амплитудно-частотная характериситка или АЧХ, а \(Arg(H_\Phi(\exp (i2\pi \omega))) = \phi_\Phi(\omega)\) - фазово-частотная характеристика или ФЧХ.

Пусть \(x_n = \cos(2\pi \omega n)\), тогда \((\Phi(X))_n = A_\Phi(\omega) \cos(2\pi \omega n + \phi_\Phi(\omega))\).

9. AЧХ фильтра скользящего среднего, зависимость от длины окна.

Рассмотрим данные: объем ежемесячной продукции молока на одну корову в фунтах с января 1962 по декабрь 1975. Всего 168 точек.

setwd(“E:/Study/Magistracy/10 term/Time series (pr)/First_an_1mar/”)

library(zoo)
setwd("E:/Study/Magistracy/10 term/Time series (pr)/Reports/")
data_raw <- read.csv("milk.csv")
colnames(data_raw) <- c("date", "val")
data_raw$date <- as.Date(data_raw$date, format = '%Y/%m/%d')
data_zoo <- zoo(data_raw[,2], as.Date(data_raw$date, format = '%Y/%m/%d'))
data <- ts(data_zoo)

plot(data_zoo, main = "Верменной ряд")

Применим к нему фильтр скользящего среднего с длиной окна 12 и 24. Выбираем длину окна кратную периоду для того, чтобы убарть сезонность.

Функция сокользящего среднего

movingAverage <- function(x, n=1, centered=TRUE) {
      
      if (centered) {
            before <- floor  ((n-1)/2)
            after  <- ceiling((n-1)/2)
      } else {
            before <- n-1
            after  <- 0
      }

      s     <- rep(0, length(x))
      count <- rep(0, length(x))
      
      new <- x
      count <- count + !is.na(new)
      new[is.na(new)] <- 0
      s <- s + new

      i <- 1
      while (i <= before) {
            new   <- c(rep(NA, i), x[1:(length(x)-i)])
            
            count <- count + !is.na(new)
            new[is.na(new)] <- 0
            s <- s + new
            
            i <- i+1
      }
      
      i <- 1
      while (i <= after) {
            new   <- c(x[(i+1):length(x)], rep(NA, i))
            
            count <- count + !is.na(new)
            new[is.na(new)] <- 0
            s <- s + new
            
            i <- i+1
      }
      s/count
}
plot(data, main = "Moving average")
ma12 <- movingAverage(data_raw[,2], 12, TRUE)
lines(ma12, col = "red")

ma24 <- movingAverage(data_raw[,2], 24, TRUE)
lines(ma24, col = "blue")
legend(0, 980, legend = c("M = 12", "M = 24"), col = c("red", "blue"), lty = 1)

Функция АЧХ фильтра

afc <- function(filter, omega)
{
      k <- seq_along(filter) - 1
      h <- function(o) sum(rev(filter)*exp(-k*1i*o))
      abs(sapply(omega, h))
}
h <- 12
freq <- seq(0, pi, 0.001)
filt <- rep(1/h, h)
omega <- freq/2/pi

plot(afc(filt, freq)~ omega, type = "l", xlab = "Frequency", ylab = "Frequency responnse", xaxt = "n")
title(main = "АЧХ фильтра скользящего среднего для окна = 12")
axis(1,at = seq(0,0.5,by = 1/12),las = 2,labels=sprintf("%s",str))

По АЧХ видно, что частоты, соответствующие сезонности (1/12, 2/12 и т.д.) фильтр обращает в 0.

h= 24
filt <- rep(1/h, h)

plot(afc(filt, freq)~ omega, type = "l", xlab = "Frequency", ylab = "Frequency responnse")
title(main = "АЧХ фильтра скользящего среднего для окна = 24")

Чем больше окно, тем сильнее фильтр подавляет высокие частоты.

10. АЧХ фильтра перехода к разностям (дифференцирования)

Рассмотрим применение фильтра перехода к разностям к нашим данным. Такой фильтр увеличивает высокие частоты и гасит низкие.

h.diff <- c(1,-1)
y.diff <- stats::filter(data,h.diff,sides = 1,method = "convolution")


plot(afc(h.diff, freq)~ omega, type = "l", xlab = "Frequency", ylab = "Frequency responnse")
title(main = "АЧХ фильтра перехода к разностям")

plot(data, ylim = c(-60, 1000), main = "Применение фильтра перехода к разностям")
lines(y.diff, col = "red")
legend("topleft", legend = c("Изначальные данные", "После фильтра"), col = c("black", "red"), lty = 1:2, cex = 0.8)

11. Что такое запаздывание и отчего оно может возникать (на примере скользящего среднего)?

Рассмотрим причинный фильтр скользящего среднего, который “смотрит только в прошлое”.

plot(data, main = "Moving average")
ma12 <- movingAverage(data_raw[,2], 12, centered = TRUE)
lines(ma12, col = "red")

ma24 <- movingAverage(data_raw[,2], 12, centered = FALSE)
lines(ma24, col = "blue")
legend(0, 980, legend = c("Обычный", "Причинный"), col = c("red", "blue"), lty = 1)

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

12. Смещение при сглаживании фильтром скользящего среднего. Роль второй производной.

Рассмотрим линейный аналог фильтра \[y(a) = \int_{-\delta}^{\delta} f(a+x) w(x)dx, \] где \(\int_{-\delta}^{\delta} w(x)dx = 1\) и \(\int_{-\delta}^{\delta}x w(x)dx = 0\). ПУсть \(f\) - гладкая функция. Применим к ней фильтр и разложим результат в ряд Тейлора \[y(a) \approx f(a) \times 1 + f'(a) \times 0 + \frac{f''(a)}{2}\int_{-\delta}^{\delta}x^2 w(x)dx.\] Последнее слагаемое в формуле является смщением при приминении фильтра. Если положить \(w(x) = \frac{1}{2\delta}\), то мы получим непрерывный аналог фильтра скользящего среднего. Смещение в этом случае будет \(\frac{\delta^2 f''(x)}{6}\),то есть будет зависить от второй производной.

13. Фильтр для подавления шума. Роль нормы коэффициентов фильтра.

Рассмотрим последовательность \(x_1, \ldots, x_n\). Применим к ней фильтр \(\sum_i a_i x_i = y\). Пусть теперь последовательность с шумом \(x_1 + \epsilon_1, \ldots, x_n +\epsilon_n\), где \(D(\epsilon_i) = \sigma^2\) Применим к ней фильтр \(\sum_i a_i (x_i + \epsilon_i) = \hat{y}\). Посчитаем дисперсию оценки. \(D(\hat y) = ||A||^2 \sigma^2\), где \(A\) - вектор из коэффициентов. Норма вектора из коэффициентов говорит о том, насколько хорошо фильтра подавляет дисперсию шума. Для скользящего среднее эта норма равна \(1/n\), то есть такой фильтр достаточно хорошо подавляет шум.

14. Как связаны периодограммы ряда до применения фильтра и после применения фильтра?

Возьмем наши данные, но сделаем их более зашумленными

Периодограмма ряда до применения фильтра

n.data <- data + 70*rnorm(168)
plot(n.data)

sp <- spec.pgram(n.data,detrend =TRUE, log = "no",fast = FALSE, pad = FALSE,taper = 0, plot = FALSE)
{plot(sp,log = "no",xaxt = "n")
      axis(1,at = seq(0,0.5,by = 1/12),las = 2,labels=sprintf("%s",str))}

На преиодограмме видно, что в данных имеется большое количество шума.

Теперь применим скользящее среднее с окном 3, чтобы убрать шум и выделить периодичность.

plot(n.data, main = "Moving average")
ma <- movingAverage(n.data, 3, TRUE)
lines(ma, col = "red")

spec.pgram(ma,detrend = TRUE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

После применения фильтра периодограмма выглядит более ровной.

15. Модели данных – аддитивная и мультипликативная

Аддитивная модель: \[X_n = T_n + S_n + R_n,\] где \(T_n\) - тренд, \(S_n\) - сезонность, \(R_n\) - шум. Для аддитивной модели предполагается постоянная амплитуда периодческой компоненты и постоянная дисперсия шума.

Мультипликативная модель: \[X_n = T_n (1 + S_n)(1 + R_n).\] Амплитуда периодической компоненты растет пропорционально тренду. После логарифимирования из мультипликативной модели получим аддитивную.

16. Методы стабилизации дисперсии в разных моделях

Рассмотрим временной ряд \(X\).

Преодразования для стабилизации дисперсии шума:

  1. \(\ln(X)\). Применимо к мультипликативной модели
  2. \(\sqrt{(X)}\). Если шум пуассоновский
  3. \(\arcsin(X)\)

17. Выделение тренда у ряда с сезонностью (выбор длины окна в скользящем среднем)

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

plot(data, main = "Скользящее среднее с окном = 12")
ma12 <- movingAverage(data_raw[,2], 12, TRUE)
lines(ma12, col = "red")

18. Переход к разностям – плюсы и минусы

Плюсы:

  1. Если ряд имеет вид \(x_t=c +y_t\), где \(c\) – тренд, \(y_t\) – стационарный процесс, то при переходе к разностям получим стационарный ряд.

  2. Если переход к разностям первого порядка, то уберем линейный тренд, если второго порядка - квадратичный и т.д.

Минусы:

  1. Увеличиваем вклад шума.

19. Скользящее среднее и скользящая медиана

Если в данных имеются выбросы, то используют скользящую медиану (т.к. скользящее среднее неустойчиво к выбросам). Однако у скользящей медианы негладкий вид и она не является оценкой математического ожидания, если распределение несимметричное.

library(stats)

plot(data)
ma12 <- movingAverage(data_raw[,2], 12, TRUE)
lines(ma12, col = "red")

mm13 <- runmed(data, 13)
lines(mm13, col = "blue")
legend("topleft", legend = c("Среднее", "Медиана"), col = c("red", "blue"), lty = 1)

20. Растекание частоты в периодограмме. Подправка длины ряда для ее устранения

Растекание частоты происходит, когда частота не попадает в решетку периодограммы. Смоделируем ряд из 50 точек состоящий из косинуса с частотой 5/13.

x <- 1:50
y <- cos(2*pi*(5/13)*x)

spec.pgram(y,detrend = FALSE, log = "no",fast = FALSE, pad = TRUE, taper = 0)

Добавим теперь два нуля, чтобы длина ряда была кратна периоду.

y <- c(y, c(0,0))

spec.pgram(y,detrend = FALSE, log = "no",fast = FALSE, pad = TRUE, taper = 0)

21. Выделение тренда с помощью параметрической регрессии

model3 <- lm(data_raw[,2] ~ poly(data_raw[,1],3))
trend_lm <- model3$fitted.values

plot(data)
lines(trend_lm, col = "red")

22. Выделение тренда с помощью метода LOESS

LOESS метод основан на сглаживании с помощью построения локальных взвешенных линейных регрессий.

Рассмотрим применение алгоритмав точке t

  1. Выбираем окно h

  2. Для точек-соседеей из \(t\pm h\) строим взвешенную регрессию. Веса могут быть, например, кубические \(w_i = (1-\Delta^3)^3,\) где \(\Delta_i = \frac{t_i-t}{h}\).

data.loess <- data_raw
data.loess$index <- 1:168

fit.loess <- loess(val~ index, data = data.loess, span = 0.4, degree = 1)      
trend.loess <- predict(fit.loess)

plot(data)
lines(trend.loess, col = "red")

25. Метод разложения Classical seasonal decomposition

Зададим период T = 12

  1. Применяем скользящее среднее с длиной окна T. Получим ряд \(X_n = \widetilde{T}_n + \widetilde{S_n + N_n}\)

  2. Сглаживаем \(\widetilde{S_n + N_n}\) с маленьким окном. Получаем \(\widetilde{S_n + N_n} = \widetilde{\widetilde{S}}_n + \widetilde{\widetilde{N}}_n\).

  3. Применяем скользящее среднее с небольшим окном к ряду \(\widetilde{T}_n + \widetilde{\widetilde{N}}_n\), получаем \(\widetilde{\widetilde{T}}_n + \widetilde{\widetilde{\widetilde{N}}}_n\)

Получаем ряд \(x_n = \widetilde{\widetilde{T_n}} + \widetilde{\widetilde{S_n}} + \widetilde{\widetilde{\widetilde{N_n}}}\).

Метод seasonal decomposition представлен в R функцией decompose.

dataForStl <- ts(data, deltat = 1/12)

fit.dec <- decompose(dataForStl, type = "additive")
plot(fit.dec)

26. Метод разложения STL

Описание алгоритма взято из этой статьи http://www.wessa.net/download/stl.pdf

Хотим представить ряд в следующем виде \(Y_n = T_n + S_n + N_n.\)

Пусть у нас уже оцененные \(T^{(k)}\), \(S^{(k)}\) на k-ом шаге алгоритма.

Inner loop

  1. Detrending. \(\widetilde{Y} = Y - T^{(k)}\)

  2. Cycle-subseries smoothing

Каждый период длины \(n_p\) сглаживается \(LOESS(q), q=n_s,d=1\). Получили ряд \(C^{(k +1)}_n,\) где \(n\) меняется от \(-n_p + 1\) до \(N + n_p\).

  1. Low-Pass Filtering of Smoothed Cycle-Subseries

К \(C^{(k +1)}_\nu\) последовательно применяем фильтр сокльзящего серднего с длиной окна \(n_p\), затем еще один фильтр с дилной окна \(n_p\). Далее применяем фильтр скользящего среднего с окном 3 и в конце LOESS c параметрами d = 1, q = \(n_l\). Полчили ряд \(L^{(k+1)}_n\).

  1. Detrending of Smoothed Cycle-Subseries

Вычисляем \(S^{(k+1)}_n = C^{(k +1)}_n - L^{(k+1)}_n\)

  1. Deseasonalizing

Вычисляем \(Y_n - S_n^{(k+1)}\)

  1. Trend Smoothing

Применяем LOESS c \(q=n_t, d=1\) к ряду без сезонности, вычисленному на предыдущем шаге.

Получаем оценку тренда \(T_n^{(k+1)}\)

Для рядов с аутлаерами есть еще дополнительный outer loop.

Таким образов STL имеет следующие параметры:

  • \(n_i\) (inner) – число итераций внутреннего цикла (р).
  • \(n_o\) (outer) – число итераций внешнего цикла. У нас в данных нет аутлаеров, поэтому outer = 0
  • \(n_l\) (l.window) – сглаживающий параметр для low-pass фильтра. Возьмем равным 13. Рекомендуют брать ближайший нечетный к периоду.
  • \(n_t\) (t.window) – сглаживающий параметр для тренда. Рекомендуют брать ближайшим нечетным к (1.5*period) / (1-(1.5/s.window)). Равен 23.
  • \(n_s\) (s.window) – сглаживающий параметр сезонности. Берем 13, так как годовая периодичность и параметр должен быть нечетным
fit.stl <- stl(dataForStl, s.window = 13, l.window = 13, outer = 0, inner = 1, t.window = 23)
plot(fit.stl)

24. Оценивание поведения дисперсии шума с помощью выделения тренда

Возьмем остатки после stl

noise <- fit.stl$time.series[,3]
plot(noise)

Мы получили шум с некоторой дисперсией \(\zeta_n = \sigma(n)\varepsilon_n\), где \()\varepsilon_n\) - белый шум. Рассмотрим квадрат \(\zeta_n^2 = \sigma(n)^2\varepsilon_n^2\) и выделим тренд \(E\zeta_n^2 = \sigma(n)^2\).

plot(noise^2, type= "l")
ma.noise <- movingAverage(noise^2, 12)
lines(ma.noise, col = "red")
legend("topleft", c("noise^2", "sigma^2"), lty=c(1,1),col=c("black","red"))

23. Нахождение огибающей периодического ряда с помощью выделения тренда

plot(noise, type = "l")
lines(sqrt(ma.noise), col = "red")
lines(-sqrt(ma.noise), col = "red")