Библиотеки и вспомогательные функции

library(ggplot2)
library(ggpubr)
## Loading required package: magrittr
library(cowplot)
## 
## ********************************************************
## Note: As of version 1.0.0, cowplot does not change the
##   default ggplot2 theme anymore. To recover the previous
##   behavior, execute:
##   theme_set(theme_cowplot())
## ********************************************************
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
library(Rssa)
## Warning: package 'Rssa' was built under R version 3.6.3
## Loading required package: svd
## Loading required package: forecast
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:ggpubr':
## 
##     gghistogram
## 
## Attaching package: 'Rssa'
## The following object is masked from 'package:stats':
## 
##     decompose
library(tseries)
## Warning: package 'tseries' was built under R version 3.6.2
library(forecast)
library(urca)
library(TTR)

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
## 
##     extract
library(gdata)
## gdata: Unable to locate valid perl interpreter
## gdata: 
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location
## gdata: of a valid perl intrpreter.
## gdata: 
## gdata: (To avoid display of this message in the future, please
## gdata: ensure perl is installed and available on the executable
## gdata: search path.)
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.
## 
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.
## 
## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.
## 
## Attaching package: 'gdata'
## The following objects are masked from 'package:dplyr':
## 
##     combine, first, last
## The following object is masked from 'package:stats':
## 
##     nobs
## The following object is masked from 'package:utils':
## 
##     object.size
## The following object is masked from 'package:base':
## 
##     startsWith
library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:gdata':
## 
##     first, last
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
library(Matrix)
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
tser.align.center <- function(tser, len, def.val = NA) {
  len.diff <- len - length(tser)
  if (len.diff <= 0) {
    return(tser)
  }
  
  l.al <- len.diff %/% 2
  r.al <- len.diff - l.al
  c(rep(def.val, l.al), tser, rep(def.val, r.al))
}

tser.align.right <- function(tser, len, def.val = NA) {
  len.diff <- len - length(tser)
  if (len.diff <= 0) {
    return(tser)
  }
  
  c(rep(def.val, len.diff), tser)
}

tser.align.left <- function(tser, len, def.val = NA) {
  len.diff <- len - length(tser)
  if (len.diff <= 0) {
    return(tser)
  }
  
  c(tser, rep(def.val, len.diff))
}

tser.align <- function(tser, len, def.val = NA, type='center') {
  tser.len <- length(tser)
  len.diff <- len - tser.len
  if (len.diff <= 0) {
    return(tser)
  }
  
  switch(type, 'center'=tser.align.center, 'left'=tser.align.left, 'right'=tser.align.right)(tser, len, def.val)
}

tser.df.labelled.h <- function(lab, tser) {
  tser.len <- length(tser)
  as.data.frame(list(1:tser.len, as.numeric(tser), rep(lab, tser.len)), col.names = c("Index", "Values", "Labels"))
}

# for aligned series
tser.df.labelled <- function(tser.labs, ...) {
  Reduce(rbind, mapply(tser.df.labelled.h, tser.labs, list(...), SIMPLIFY = F, USE.NAMES = F))
}

tser.my.plot.h <- function(tser.df) {
  ggplot(data = tser.df) + geom_line(aes(x = Index, y = Values, col = Labels), alpha = 0.75)
}

tser.my.plot <- function(tser.orig, tser.trend, tser.season=NULL, plot.resid=T, plot.season=T, plot.reco=T,
                         component.labs=c("Trend", "Seasonal"), type.align='center', align.def.val=NA) {
  tser.trend.al <- tser.align(tser.trend, length(tser.orig))
  tser.resid <- tser.orig - tser.trend.al
  
  plot.list <- list(tser.my.plot.h(tser.df.labelled(c("Original", component.labs[1]), tser.orig, tser.trend.al)))

  if (!is.null(tser.season) && plot.season) {
    tser.season.al <- tser.align(tser.season, length(tser.orig))
    plot.list[[length(plot.list) + 1]] <- tser.my.plot.h(tser.df.labelled(c("Detrend", component.labs[2]), tser.resid, tser.season.al))
    tser.resid <- tser.resid - tser.season.al
  }
  
  if (!is.null(tser.season) && plot.reco) {
    tser.season.al <- tser.align(tser.season, length(tser.orig))
    plot.list[[length(plot.list) + 1]] <-  tser.my.plot.h(tser.df.labelled(c("Original", "Reconstructed"), tser.orig,
                                                                           tser.trend.al + tser.season.al))  
  }
  
  if (plot.resid) {
    plot.list[[length(plot.list) + 1]]<- tser.my.plot.h(tser.df.labelled(c("Residual"), tser.resid))
  }
  
  ggarrange(plotlist = plot.list, ncol = 1, nrow = length(plot.list), align = 'v')
}
train.test.split <- function(time.series, test.size, revert=F) {
  if(revert) {
      time.series <- replace(time.series, TRUE, rev(time.series))
  }
  list(Train = head(time.series, length(time.series) - test.size), Test = tail(time.series, test.size))
}

Данные для предсказания

В этом отчете мы будем рассматривать способы предсказания будущих значений. Как прежде, мы продолжим рассматривать данные о доходах компании Johnson & Johnson:

plot(JohnsonJohnson)

Чтобы оценить качество прогноза, мы разделим ряд на две части: по первой мы будем строить модель, а со второй частью мы будем сравнивать полученные предсказания. Давайте мы возьмем для построения модели первые \(80\) наблюдений, а для сравнения – последние \(4\). Так мы постараемся предсказать доход компании на следующий год

test.size <- 4
end.split <- train.test.split(JohnsonJohnson, test.size)

Прогнозирование с помощью экспоненциального сглаживания (ETS)

Аббревиатура \(ETS(\cdot, \cdot, \cdot)\) расшифровывается как Error, Trend, Seasonal. Далее мы расскажем, что модель может быть применена для ряда случаев: для аддитивных и мультипликативных ошибок, постоянного или монотонного тренда, мультипликативной или аддитивной сезонности. Каждый из случаев носит название, состоящее из трех букв. Например ETS(M, A, N) будет означать, что у ряда мультипликативные ошибки и монотонный тренд без сезонных колебаний.

Простое экспоненциальное сглаживание

В первом отчете мы рассматривали скользящее среднее. Мы можем использовать его для построения простейшего предсказания: \(\hat x_{N+1} = \frac{1}{N} \sum_{i=1}^{N}{x_i}\), где \(N\) – длина исходного ряда. Заметим, что в таком определении все точки получают одинаковый вес, одинаково влияют на предсказание. Нам бы хотелось, чтобы влияние старых наблюдений было меньше.

Давайте сделаем такой фильтр. Выберем некоторое \(0 < \alpha < 1\) и положим \(\hat x_{N+1} = \alpha x_{N} + \alpha (1 - \alpha) x_{N-1} + \alpha (1 - \alpha) ^ 2 x_{N-2} + \ldots\). При таком определении у старых наблюдений вес убывает экспоненциально.

Перепишем нашу формулу: \[\hat x_{N+1} = \alpha x_N + (1 - \alpha) l_{N}\] Аналогично мы определим и предыдущие элементы ряда: \[2 \leq n \leq N,\ l_{n} = \alpha x_{n-1} + (1 - \alpha) l_{n-1}\]

На самом деле мы определили бесконечный линейный фильтр (Infinite Response Filter), но мы имеем дело с конечными последовательностями, поэтому нам необходимо найти начальное значение, дальше которого формулу мы не разворачиваем: \(l_2 = \alpha x_1 + (1 - \alpha) l_1\). \(l_1\) – это начальное значение для нашего фильтра, и мы находим его с помощью метода наименьших квадратов вместе с параметром \(\alpha\). Теперь мы можем явно выразить предсказание: \[\hat x_{N+1} = \sum_{i=1}^{N}{\alpha (1-\alpha) ^ {N-i} x_{i}} + (1 - \alpha) l_1\]

Для дальнейших рассуждений нам удобно представить формулу в составном виде: \[\begin{align*}\hat x_{N+t} &= l_{N+1}, \ t \geq 1 \\ l_{n+1} &= \alpha x_n + (1 - \alpha) l_n, 1 \leq n \leq N \end{align*}\] \(x_n\) – данные, с которыми мы работаем, исходный временной ряд. \(l_n\) – предсказанное для \(x_n\) значение с помощью сглаживания.

С помощью такой формулы мы сможем предсказывать значение только для рядов без тренда и сезонности: мы предсказываем одно и то же будущее значение и никак не учитываем изменения ряда.

Выделение линейного тренда

Хольтом (Holt) в 1957 было предложено расширение метода для рядов с линейным трендом:

\[\begin{align*}\hat x_{N+t} &= l_{N+1} + t \cdot b_{N+1},\ t \geq 1 \\ l_{n+1} &= \alpha x_n + (1 - \alpha)(l_n - b_n), 1 \leq n \leq N \\ b_{n+1} &= \beta^*(l_{n+1} - l_n) + (1 - \beta^*) b_n, 1 \leq n \leq N \end{align*}\]

\(b_{N+1}\) определяет скорость роста ряда, также как и для \(l_n\), нам понадобится найти начальное значение параметра \(b_1\) и \(\beta^*\) с помощью метода наименьших квадратов.

Замедление скорости роста

Описанная ранее формула предполагает неограниченный рост значений ряда для \(t \rightarrow \infty\). Мы хотели бы предсказывать не столь большие значения для далеких предсказаний. Для этого мы немного изменим формулу:

\[\begin{align*} \hat x_{N+t} &= l_{N+t} + (\sum_{i=1}^{t}{\phi^t}) b_{N+1} \\ l_{n+1} &= \alpha x_n + (1 - \alpha)(l_n - \phi b_n), 1 \leq n \leq N \\ b_{n+1} &= \beta^*(l_{n+1} - l_n) + (1 - \beta^*) \phi b_n, 1 \leq n \leq N \end{align*}\]

Полагая \(0 < \phi < 1\), получим что сумма \(\sum_{i=1}^{t}{\phi^t}\) оказывается конечной для \(t \rightarrow \infty\) и равна \(\frac{1}{1 - \phi}\). Таким образом мы можем ограничить значения для далеких предсказаний.

Holt-Winters сезонное сглаживание

Винтерс (Winters) в 1960 расширил метод, предложенный Хольтом, добавив сезонную составляющую. Обозначим \(m\) – частоту сезонности, \(m \geq 1\). Например, для ежемесячных наблюдений \(m=12\), для квартальных – \(m=4\).

\[\begin{align*}\hat x_{N+t} &= l_{N+1} + t \cdot b_{N+1} + s_{N + t - m (\lfloor \frac{t-1}{m} \rfloor + 1)},\ t \geq 1 \\ l_{n+1} &= \alpha (x_n - s_{n - m + 1}) + (1 - \alpha)(l_n - b_n), 1 \leq n \leq N \\ b_{n+1} &= \beta^*(l_{n+1} - l_n) + (1 - \beta^*) b_n, 1 \leq n \leq N \\ s_{n+1} &= \gamma(x_{n+1} - l_n - b_n) + (1 - \gamma)s_{n - m + 1} \end{align*}\]

Добавилось ещё \(m + 1\) параметров: \(s_1, \ldots, s_m\) – начальные значения для сезонностей и параметр \(\gamma\), которые нам нужно определить как и раньше с помощью МНК. Мы также можем замедлить скорость роста тренда при \(t \rightarrow \infty\), заменив \(t \cdot b_{N+1}\) на \(\sum_{i=1}^{t}{\phi^t} b_{N+1}\).

Также метод может быть обобщен для мультипликативных рядов.

State Space Model

До сих пор мы рассматривали методы для точечных предсказаний, но мы хотим также строить доверительные интервалы для них, а также уметь сравнивать разные методы. Нам понадобится ввести статистическую модель для этого.

Для начала рассмотрим простое экспоненциальное сглаживание:

\[\begin{align*}\hat x_{N+1} &= l_{N+1} \\ l_{n+1} &= \alpha x_n + (1 - \alpha) l_n \end{align*}\]

Мы можем переписать второе уравнение в следующем виде: \[l_{n+1} = l_n + \alpha(x_n - l_n)\]

Заметим, что второе слагаемое в правой части соответствует ошибке предсказания, давайте обозначим её за \(\epsilon_n\). Тогда \(l_{n+1} = l_n + \alpha \epsilon_n\). Предположим, что \(\epsilon_n \sim N(0, \sigma^2)\), независимы и одинаково распределенные.

Мы получим следующие уравнения:

\[\begin{align*} x_n &= l_n + \epsilon_n \\ l_{n+1} &= l_n + \alpha\epsilon_n \end{align*}\]

Первое позволяет установить связь с элементами исходного ряда, а последнее задаёт правило перехода к предсказанию следующего наблюдения. Заметим, что если \(\alpha = 1\), то мы получим простое случайное блуждание.

Для модели с линейным трендом мы обозначим \(\epsilon_n = x_n - l_n - b_n\), тогда мы можем переписать наши формулы в следующем виде:

\[\begin{align*} x_n &= l_n + b_n + \epsilon_n \\ l_{n+1} &= l_n + b_n + \alpha\epsilon_n \\ b_{n+1} &= b_n + \beta \epsilon_n \end{align*}\] Где \(\beta = \alpha \beta^*\). Для ряда с сезонностью с периодом \(m\) и линейным трендом, предполагая \(\epsilon_n = x_n - l_n - b_n - s_{n - m}\), мы получим следующие формулы:

\[\begin{align*} x_n &= l_n + b_n + s_{n - m + 1} + \epsilon_n \\ l_{n+1} &= l_n + b_n + \alpha\epsilon_n \\ b_{n+1} &= b_n + \beta \epsilon_n \\ s_{n+1} &= s_{n-m+1} + \gamma \epsilon_n \end{align*}\]

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

Выбор модели

Теперь, когда мы определили модель данных, мы можем находить параметры с помощью метода максимального правдоподобия. Кроме того, с помощью метода максимального правдоподобия мы можем сравнивать разные модели.

Определим информационные критерии на основе значения функции правдоподобия:

\[\begin{align*} AIC &= -2log(L) + 2k \\ AIC_c &= AIC + \frac{k(k+1)}{N - k -1} \\ BIC &= AIC + k[log(N) - 2]\end{align*}\]

Где \(L\) – значение функции правдоподобия, \(k\) – количество параметров модели. Такие критерии позволяют учитывать сложность модели и количество наблюдений.

Применение к реальным данным

Предсказание для конца ряда

Мы хотим определить модель для нашего ряда:

plot(end.split$Train)

Среди моделей ETS есть те, которые предназначены для рядов с мультипликативной сезонностью и мультипликативными ошибками, поэтому нам не нужно пытаться сводить модель к аддитивной. Для того, чтобы определить параметры модели, мы воспользуемся функцией ets из пакета forecast. Если не ограничивать модели (например с помощью параметра additive.only), то будут рассмотрены все возможные модели и выбраны с наибольшим значением информационного критерия.

jj.ets <- ets(end.split$Train, ic = "aicc", opt.crit = "lik")
jj.ets
## ETS(M,A,A) 
## 
## Call:
##  ets(y = end.split$Train, opt.crit = "lik", ic = "aicc") 
## 
##   Smoothing parameters:
##     alpha = 0.1114 
##     beta  = 0.1111 
##     gamma = 0.3604 
## 
##   Initial states:
##     l = 0.6246 
##     b = -0.0051 
##     s = -0.1659 0.165 -0.0087 0.0095
## 
##   sigma:  0.0932
## 
##      AIC     AICc      BIC 
## 142.2030 144.7744 163.6412

Модель ETS(M, A, A) с мультипликативными остатками, аддитивным трендом и сезонностью.

Сделаем предсказание:

jj.ets.forecast <- forecast(jj.ets, h = test.size)
plot(jj.ets.forecast)
lines(JohnsonJohnson)

Предсказание повторяет форму ряда. Давайте сравним с исходным рядом:

tser.my.plot(end.split$Test, jj.ets.forecast$mean, component.labs = c("Predicted"))

Предсказание повторяет форму ряда. Остатки заметно отличаются от нуля. Проверять, соответствуют ли они предположению об i.i.d мы не будем.

Среднее квадратов отклонений:

jj.ets.mse <- mean((end.split$Test - jj.ets.forecast$mean) ^ 2)
jj.ets.mse
## [1] 0.4987922

Сейчас мы не использовали подавление тренда для предсказания. Мы делаем прогноз всего на четыре значения вперед и, по идее, подавление не должно оказывать большого влияния, но все же проверим:

jj.ets.d <- ets(end.split$Train, ic = "aicc", opt.crit = "lik", damped = T)
jj.ets.d
## ETS(M,Ad,A) 
## 
## Call:
##  ets(y = end.split$Train, damped = T, opt.crit = "lik", ic = "aicc") 
## 
##   Smoothing parameters:
##     alpha = 0.084 
##     beta  = 0.0817 
##     gamma = 0.6881 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 0.5691 
##     b = 0.0376 
##     s = -0.2754 0.2574 -0.0337 0.0516
## 
##   sigma:  0.1
## 
##      AIC     AICc      BIC 
## 153.3716 156.5600 177.1918

Заметим, что демпирующий параметр \(\phi = 0.98\), то есть подавление имеет значение. Для самого последнего предсказания мы получим коэффициент для тренда: \(\phi + \phi^2 + \phi^3 + \phi^4 \approx 3.8\), вместо \(4\) для модели без демпирования. Сделаем предсказание:

jj.ets.d.forecast <- forecast(jj.ets.d, h = test.size)
plot(jj.ets.d.forecast)
lines(JohnsonJohnson)

Кажется, что вторая точка лежит чуть ниже. Давайте сравним с исходным рядом:

tser.my.plot(end.split$Test, jj.ets.d.forecast$mean, component.labs = c("Predicted"))

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

jj.ets.d.mse <- mean((end.split$Test - jj.ets.d.forecast$mean))
jj.ets.d.mse
## [1] 0.1765835

Действительно, MSE заметно меньше, чем для модели без демпирования.

Прогнозирование с помощью ARIMA

Ряд, у которого свойства не зависят от времени, называется стационарным. Любой ряд с трендом или сезонностью стационарным не является.

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

Нам удобно ввести оператор сдвига \(B\): \(Bx_n = x_{n-1}\). Мы покажем, как с помощью такого оператора в наглядно форме представить модель ARIMA.

Модель AR(p)

Авторегрессионной моделью степени \(p\) мы назовем случайный процесс, значения которого задаются по формуле: \[x_n = c + \sum_{i=1}^{p}{\phi_i x_{n-i}} + \epsilon_n\]

Где \(\epsilon_n \sim N(0, \sigma^2)\) – независимые случайные величины.

Рассмотрим модели, возможные для AR(1):

  • Если \(\phi_1 = 0,\ c=0\) – то модель соответствует белому шуму
  • Если \(\phi_1 = 1,\ c=0\) – то модель соответствуем случайному блужданию
  • Если \(\phi_1 = 1,\ c \neq 0\) – то модель соответствуем случайному блужданию с началом в точке \(c\)
  • Если \(\phi < 0\), то значения \(x_n\) будут распределены вокруг среднего значения

Заметим, что \(x_n\) коррелирует со всеми \(x_i\) для \(i < n\): последовательно заменяя \(\forall i < n: x_{i}\) по формулам, мы получим зависимость \(x_i\) от \(x_i\) для \(i < n - p\). На самом деле, корреляции будут убывать с увеличением лага.

С другой стороны, только первые \(p\) частных автокорреляций будут отличны от нуля (Это можно объяснить тем, что явно \(x_n\) зависит только от предыдущих \(p\) значений, причем линейно. Логично предположить, что корреляция \(x_n\) с \(x_{n-p-1}\), за вычетом влияния \(x_{n-1},\ldots,x_{n-p}\), должна быть неотличима значимо от нуля)

Мы можем записать модель с помощью оператора сдивга: \((1 - \phi_1 B - \phi_2B^2 - \ldots - \phi_pB^p)x_n = c + \epsilon_n\).

Модель MA(q)

Моделью скользящего среднего степени \(q\) мы назовем модель \[x_n = c + \epsilon_n + \sum_{i=1}^{q}{\theta_i \epsilon_{n-i}}\] Здесь \(\epsilon_n \sim N(0, 1)\), независимые и одинаково распределенные случайные величины.

Для процесса, который описывается такой моделью, \(x_n\) будет коррелировать только с предыдущими \(q\) эелементами в силу независимости \(\epsilon_n\): как только \(x_i\) перестаёт включать \(\epsilon_n, \epsilon_{n-1}, \ldots, \epsilon_{n-q}\), \(\rho(x_i, x_n) \approx 0\). Только первые \(q\) автокорреляций будут значимо отличаться от нуля. С другой стороны, частные автокорреляции будут убывать постепенно.

Эту модель мы также можем записать с помощью оператора сдвига: \(x_n = c + (1 + \theta_1 B + \theta_2 B^2 + \ldots + \theta_q B^q )\epsilon_n\).

Модель ARMA(p, q)

Мы можем объединить предыдущие модели, чтобы описать более широкий класс стационарных процессов. Легко выразить такую модель с помощью оператора сдвига: \[(1 - \phi_1 B - \phi_2 b^2 - \ldots - \phi_p B^p) x_n = c + (1 + \theta_1 B + \theta_2 B^2 + \ldots + \theta_q B^q) \epsilon_t\]

В левой части стоит AR(p), а в правой – MA(q).

Нахождение параметров модели

В случае, если мы имеем дело с чистой AR(p) или MA(q) моделью, то могут помочь графики частных или обычных автокорреляций соответственно. В таком случае, номер последнего значимо отличающегося от нуля коэффициента скорее всего и будет порядком модели.

Если мы имеем дело со смешанной моделью ARMA(p, q), то графики уже не помогут. Для того, чтобы определить \(p\) и \(q\) мы можем воспользоваться информационными критериями. Обозначим \(\Theta = (\phi_1, \phi_2, \ldots, \phi_p, \theta_1, \theta_2, \ldots, \theta_q)\) – параметры модели, \(X_N\) – ряд, \(L = L(\Theta | X_N)\) – максимум функции правдоподобия для модели со степенями \(p\) и \(q\). Мы рассмотрим несколько информационных критериев:

\[\begin{align*} AIC &= -2 \log{(L)} + 2(p + q + k + 1) \\ AIC_c &= AIC + \frac{2(p + q + k + 1)(p + q + k + 2)}{N - p - q - k - 2} \\ BIC &= AIC + \lfloor \log{(N)} - 2\rfloor(p + q + k + 1) \end{align*}\]

Мы выбираем какой-нибудь из критериев, перебираем несколько пар значений (p, q) и останавливаемся на той модели, для которой значение выбранного критерия наименьшее.

Сведение процесса к стационарному

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

Пусть \(x_N\) – исходный ряд, тогда \(x^\prime_{N-1} = (x_2 - x_1, x_3 - x_2, \ldots, x_N - x_{N-1})\) – ряд из конечных разностей первого поядка. \(x^{\prime\prime}_{N-2} = (x^\prime_2 - x^\prime_1, x^\prime_3 - x^\prime_2, \ldots, x^\prime_{N-1} - x^\prime_{N-2})\) – ряд из разностей второго порядка. Операцию перехода к разностям мы назовем дифференцированием.

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

Операцию дифференцирования мы также можем выразить с помощью оператора сдвига: \(x^\prime_n = (1 - B)x_n\). Очень легко выражаются и производные высших порядков: \(x^{(d)}_n = (1 - B)^dx_n\)

ARIMA(p, d, q)

Мы разобрались с моделями для стационарных процессов. Теперь мы хотим построить модель для ряда с трендом. С помощью разностей мы сведем процесс к стационарному и будем искать модель уже для него. В итоге модель представляется так:

\[(1 - \phi_1 B - \phi_2 B^2 - \ldots - \phi_p B^p) (1 - B)^dx_n = c + (1 + \theta_1 B + \theta_2 B^2 + \ldots + \theta_q B^q) \epsilon_t\]

В левой части добавился множитель \((1-B)^d\) соответствующий дифференцированию \(d\) порядка. Отметим, что информационные критерии не подходят для определения \(d\), только для поиска \(p\) и \(q\).

Проверка стационарности

Теперь нам нужно определить значение \(d\). Если у нас есть критерий стационарности, ты мы можем дифференцировать ряд до тех пор, пока он не станет стационарным, определяя момент остановки с помощью теста. Рассмотрим правую часть модели: \(\Phi(B) = (1 - \phi_1 B - \phi_2 B^2 - \ldots - \phi_p B^p) (1 - B)^d\). Если есть корни такого полинома, равные по модулю единице, то мы имеем дело с нестационарным процессом. Мы приведем примеры тестов, которые так ли иначе связаны с проверкой единичных корней (они называются Unit-Root tests):

  • Kwiatkowski, Phillips, Schmidt, Shin test (KPSS test), \(H_0\) – ряд стационарный, \(H_1\) – есть единичный корень
  • Dickey–Fuller test, \(H_0\) – есть единичный корень (ряд не стационарен)
  • Phillips–Perron test, \(H_0\) – есть единичный корень

ARIMA для ряда с сезонностью

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

Пусть частота сезонности равна \(m\). Мы можем рассматривать \(m\), например, как количеств наблюдений за год. Тогда модель ARIMA с учетом сезонности будет выглядеть так:

\[(1 - \phi_1 B - \phi_2 B^2 - \ldots - \phi_p B^p) (1 - \Phi_1 B^m - \Phi_2 B^{2m} - \ldots - \Phi_P B^{mP}) (1 - B)^d (1 - B^m)^D x_n =\\c + (1 + \theta_1 B + \theta_2 B^2 + \ldots + \theta_q B^q) (1 + \Theta_1 B^m + \Theta_2 B^{2m} + \ldots + \Theta_Q B^{mQ}) \epsilon_t\]

Мы можем ввести сезонное дифференцирование: \((1 - B^m)x_n\). Если ряд один раз сезонно-дифференцируем, то такая операция сведет его к стационарному процессу. Если ряд имеет и сезонность, и тренд, то двойное дифференцирование \((1-B)(1-B^m)x_n\) приведет к стационарности.

Определение параметров для модельного ряда

model.ts <- read.ts('./arima/ts11.txt', header = T)
model.ts.train <- window(model.ts, start = 1, end = 900)
model.ts.test <- window(model.ts, start = 901, 1000)
plot(window(model.ts.train, 0, 100))
## Warning in window.default(x, ...): 'start' value not changed

График похож на стационарный процесс: нет явного тренда или периодичности. Давайте применим тест ADF (Augumented Dickey-Fuller), чтобы проверить гипотезу о наличии единичного корня:

adf.test(model.ts.train, alternative = "stationary")
## Warning in adf.test(model.ts.train, alternative = "stationary"): p-value
## smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  model.ts.train
## Dickey-Fuller = -9.3026, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

Гипотеза о том, что есть единичный корень, не подтвердилась, нам не нужно дифференцировать наш ряд. Давайте посмотрим на графики автокорреляций

ggAcf(model.ts.train)

ggPacf(model.ts.train)

Скорее всего, мы имеем дело со смешанной моделью. Мы не видим на графиках автокорреляций поведения, характерного для чистой AR(p) или MA(q) модели: на обоих графиках только два первых коэффициента значимо отличаются от нуля, коэффициенты убывают резко.

Попробуем подобрать модель с помощью информационных критериев. Мы не будем ускорять поиск, и просмотрим все модели (stepwise=F), мы ограничимся моделями только для стационарных процессов с нулевым средним (stationary=T, allowmean=F). Мы будем определять модель по AICc:

auto.arima(model.ts.train, stepwise = F, stationary = T, seasonal = F, allowmean = F, ic = "aicc", trace = T, max.order = 4)
## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(0,0,0) with zero mean     : 1803.56
##  ARIMA(0,0,1) with zero mean     : 1744.222
##  ARIMA(0,0,2) with zero mean     : 1735.387
##  ARIMA(0,0,3) with zero mean     : 1736.392
##  ARIMA(0,0,4) with zero mean     : 1737.838
##  ARIMA(1,0,0) with zero mean     : 1763.691
##  ARIMA(1,0,1) with zero mean     : 1739.035
##  ARIMA(1,0,2) with zero mean     : 1737.772
##  ARIMA(1,0,3) with zero mean     : 1739.203
##  ARIMA(2,0,0) with zero mean     : 1736.954
##  ARIMA(2,0,1) with zero mean     : 1736.897
##  ARIMA(2,0,2) with zero mean     : 1737.689
##  ARIMA(3,0,0) with zero mean     : 1736.66
##  ARIMA(3,0,1) with zero mean     : 1737.977
##  ARIMA(4,0,0) with zero mean     : 1738.446
## 
##  Now re-fitting the best model(s) without approximations...
## 
## 
## 
## 
##  Best model: ARIMA(0,0,2) with zero mean
## Series: model.ts.train 
## ARIMA(0,0,2) with zero mean 
## 
## Coefficients:
##          ma1      ma2
##       0.2586  -0.1109
## s.e.  0.0334   0.0334
## 
## sigma^2 estimated as 0.4009:  log likelihood=-864.72
## AIC=1735.45   AICc=1735.48   BIC=1749.86

auto.arima советует выбрать модель MA(2). Отметим, что значения \(AIC_c\) для \(ARIMA(2, 0, 0)\) и \(ARIMA(0, 0, 2)\) довольно близки. На самом деле, если взять весь ряд целиком, то auto.arima предложит модель AR(2). Дальше для предсказаний мы используем модель MA(2).

Прогноз для модельных данных

model.arima <- auto.arima(model.ts.train, stepwise = F, stationary = T, seasonal = F, allowmean = F, ic = "aicc", trace = F, max.order = 4)
model.forecast <- forecast(model.arima, h = 100)
plot(model.forecast)

Раз мы имеем дело со стационарным рядом, то и предсказания не должны меняется со временем: модель предсказывает ноль. Посмотрим на остатки:

ggtsdisplay(model.ts.test - model.forecast$mean)

Действительно, судя по графикам автокорреляций мы имеем дело с белым шумом: у нас нет значимых коэффициентов. Мы правильно определили модель.

Применение к реальным данным

Наши данные выглядят так:

plot(end.split$Train)

Периодичность наших данных \(m = 4\). Мы попробуем применить ARIMA как к исходным данным, так и после преобразования Бокса-Кокса.

Применение к исходным данным

В принципе, по графику видно, что ряд не стационарный. Подберем \(d\):

ndiffs(end.split$Train, test = "adf")
## [1] 1

И для сезонной составляющей:

nsdiffs(end.split$Train)
## [1] 1

Давайте посмотрим, как выглядит сезонно-продифференцированный ряд:

ggtsdisplay(diff(diff(end.split$Train), lag = 4))

Остатки явно гетероскедастичны. Нам бы все-таки хотелось, чтобы дисперсия была однородна. Возможно, преобразование Бокса-Кокса позволит решить эту проблему.

Подберем оставшиеся параметры ARIMA:

jj.arima <- auto.arima(end.split$Train, d = 1, D = 1, ic = "aicc", stepwise = F, trace = F)
jj.arima
## Series: end.split$Train 
## ARIMA(3,1,1)(0,1,0)[4] 
## 
## Coefficients:
##           ar1     ar2      ar3      ma1
##       -0.2213  0.1175  -0.3158  -0.5893
## s.e.   0.1917  0.1821   0.1254   0.1847
## 
## sigma^2 estimated as 0.176:  log likelihood=-39.84
## AIC=89.68   AICc=90.55   BIC=101.27

Предскажем значения:

jj.arima.forecast <- forecast(jj.arima, h = test.size)
plot(jj.arima.forecast)
lines(JohnsonJohnson)

Сравним с исходным рядом:

tser.my.plot(end.split$Test, jj.arima.forecast$mean, component.labs = c("Predicted"))

Предсказания повторяют форму ряда, остатки распределены около нуля.

Сумма квадратов ошибок:

jj.arima.mse <- mean((end.split$Test - jj.arima.forecast$mean) ^ 2)
jj.arima.mse
## [1] 0.2710524

После преобразования Бокса-Кокса

Подберем параметр \(\lambda\) преобразования Бокса-Кокса:

jj.lambda <- BoxCox.lambda(end.split$Train)
jj.lambda
## [1] 0.1803894

Посмотрим на преобразованные данные:

jj.log <- BoxCox(end.split$Train, jj.lambda)
plot(jj.log)

Разница в дисперсии наблюдений стала чуть меньше, но все же полностью не ушла. Мы видим, что для наблюдений из середины отрезка разброс явно меньше, чем на концах.

Найдём параметры \(d\) и \(D\):

ndiffs(jj.log, test = "adf")
## [1] 1
nsdiffs(jj.log)
## [1] 1

Посмотрим на то, как выглядит ряд после дифференцирования и на графики автокорреляций:

ggtsdisplay(diff(diff(jj.log), 4))

Дифференцированный ряд выглядит лучше, чем в предыдущем эксперименте: почти нет изменения дисперсии с течением времени.

Подберем оставшиеся параметры ARIMA:

jj.log.arima <- auto.arima(jj.log, d = 1, D = 1, stepwise = F, seasonal = T, ic = "aicc")
jj.log.arima
## Series: jj.log 
## ARIMA(0,1,1)(0,1,1)[4] 
## 
## Coefficients:
##           ma1     sma1
##       -0.6898  -0.3057
## s.e.   0.0957   0.1197
## 
## sigma^2 estimated as 0.01193:  log likelihood=60.06
## AIC=-114.13   AICc=-113.79   BIC=-107.17

На предлагается взять модель ARIMA(0, 1, 1)(0, 1, 1)

Сделаем предсказание:

jj.log.forecast <- forecast(jj.log.arima, h = test.size)
plot(jj.log.forecast)
lines(jj.log)

Сделаем обратное преобразование Бокса-Кокса и сравним с тестовыми данными:

jj.log.forecast.inv <- InvBoxCox(jj.log.forecast$mean, lambda = jj.lambda)
tser.my.plot(end.split$Test, jj.log.forecast.inv)

Кажется, что разница незначительна. Посмотрим на среднеквадратичную ошибку:

jj.log.arima.mse <- mean((end.split$Test - jj.log.forecast.inv) ^ 2)
jj.log.arima.mse
## [1] 0.145375

MSE оказался меньше.

Прогнозирование с помощью SSA

Пусть у нас есть сигнал \(S_N\) длины \(N\). Зафиксируем длину окна \(L\), \(\mathcal L_L (S_N)\) – линейное пространство порожденное сигналом \(S_N\) с окном \(L\). Число \(s_{N+1}\) называется точным L-продолжением, если отрезок длины \(L\) продолженного сигнала \((s_{N-L+1}, \ldots, s_{N+1}) \in \mathcal L (S_N)\).

Сигнал \(S_N\) L-продолжим, если L-продолжение существует и единственно. Из замечаний стоит отметить, что если ряд продолжим на один элемент вперед, то он продолжим на произвольное количество элементов вперед, и что условие продолжимости довольно сильное: очень редко когда сигнал L-продолжим.

Мы рассмотрим два способа продолжения ряда: рекуррентный и векторный.

Рекуррентный подход

Рекуррентный подход использует линейную рекуррентную формулу для построения продолжения.

Пусть на вход мы получили ряд \(X_N\) длины \(N\) и хотим построить прогноз на \(M\) точек. Также предположим, что \(X_N = S_N + \mathcal E_N\), и сигнал приближенно сильно отделим от шума.

Зафиксируем \(L\), построим L-траекторную матрицу нашего ряда. Пусть \(\mathcal L^r_L \subset \mathbb R ^ L\) – линейное пространство, порожденное сигналом размерности \(r < L\). \(P_1, \ldots, P_r\) – ортонормированный базис этого пространства. (Мы можем взять \(P_i\) из SVD разложения траекторной матрицы).

Предположим, что \(e_L \notin \mathcal L^r_L\) (\(e_L\) – L-ый орт пространства \(\mathbb R^L\)). Тогда мы можем построить линейную рекуррентную формулу порядка \(L-1\) с коэффициентами: \[(a_{L-1}, a_{L-2}, \ldots, a_1) = \frac{1}{1 - \nu^2} \sum_{i=1}^{r}{\pi_i \underline{P}_i}\]

Где \(P_i = (p_1, \ldots, p_{L-1}, \pi_i) ^ \intercal\), \(\nu^2 = \sum_{i}^{r}{\pi_i^2}\). Если выполняется условие \(e_L \notin \mathcal L_L^r\), то \(\nu^2 < 1\), а вектор \((a_{L-1}, a_{L-2}, \ldots, a_1, 1)\) принадлежит ортогональному дополнению пространства \(\mathcal L_L^r\).

Пусть \(\tilde{s}_1, \tilde{s}_2, \ldots, \tilde{s}_N\) – восстановленный сигнал. Мы построим продолженный ряд \(G_{N+M}\), такой что \[g_i = \begin{cases} \tilde{s}_i,& \mbox{if } i \leq N \\ \sum_{j=1}^{L-1}{a_j g_{i - j}},& \mbox{if }i > N \end{cases}\]

Векторный подход

Заметим, что продолжение ряда рекуррентным способом, вообще говоря, можем не лежать в пространстве \(\mathcal L^r_L\). Пусть как и раньше мы получаем на вход ряд \(X_N\), фиксируем размер окна \(L\), группируем тройки из SVD разложения траекторной матрицы и получаем восстановленную по тройкам матрицу сигнала \(\hat {\mathbb S} = [\hat S_1:\ldots: \hat S_K]\). Эта матрица ещё до диагонального усреднения и может не быть ганкелевой. Мы хотим сделать прогноз на \(M\) точек вперед.

Мы хотим построить ещё \(M\) векторов \(Z_{K+1}, \ldots, Z_{K+M} \in \mathbb R ^ L\), таких что \(Z_i \in colspan(\hat{\mathbb S})\), а матрица \(\hat {\mathbb G} = [\hat S_1:\ldots:\hat S_K : Z_{K+1} : \ldots : Z_{K+M}]\) – приближенно ганкелева.

Имея приближенно ганкелевую матрицу \(\hat {\mathbb G}\) с помощью диагонального мы перейдем к матрице \(\mathbb G \in \mathbb R ^ {L \times K + M}\), которая уже соответствует ряду \(G_{N+M}\). Последние \(M\) элементов этого ряда могут рассматриваться в качестве прогноза.

Применение к реальным данным

Восстановление сигнала

jj.ssa <- ssa(end.split$Train)
plot(jj.ssa)

Судя по собственным числам, нам достаточно взять первые шесть троек для сигнала (именно такое количество мы брали в предыдущих отчетах). Посмотрим на w-корреляции:

plot(jj.ssa, type = 'wcor')

plot(jj.ssa, type = 'vectors')

В принципе, седьмая компонента уже заметно смешалась с остатком, попробуем взять первые шесть для восстановления сигнала. На самом деле, шестая неустойчиво выделилась и поэтому посмотрим на прогноз по первым пяти.

R-SSA прогноз для конца ряда

jj.rssa.forecast <- forecast(jj.ssa, groups = list(Signal = c(1:5)), len = test.size,
                             method = "recurrent", interval = "confidence", only.intervals = T, level = 0.99)
plot(jj.rssa.forecast)
lines(JohnsonJohnson)

Две точки прогноза очень далеко лежат от истинных значений, причем доверительный интервал даже не накрыл истинных значений для этих точек.

Сравним только предсказания:

tser.my.plot(end.split$Test, jj.rssa.forecast$mean, component.labs = c("Predicted"))

Все-таки стоит отметить, что форма прогноза повторяет форму ряда. Давайте посчитаем среднеквадратичную ошибку:

jj.rssa.mse <- mean((end.split$Test - jj.rssa.forecast$mean) ^ 2)
jj.rssa.mse
## [1] 1.279987

Сравнительно большая ошибка.

V-SSA прогноз

jj.vssa.forecast <- forecast(jj.ssa, groups = list(Signal = 1:5), len = test.size,
                             method = "vector", interval = "prediction", only.intervals = T, level = 0.95)
plot(jj.vssa.forecast)
lines(JohnsonJohnson)

Доверительный интервал хотя бы накрывает истинные значения.

Сравним только с предсказанным отрезком:

tser.my.plot(end.split$Test, jj.vssa.forecast$mean, component.labs = c("Predicted"))

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

Посмотрим на среднеквадратичную ошибку:

jj.vssa.mse <- mean((jj.vssa.forecast$mean - end.split$Test) ^ 2)
jj.vssa.mse
## [1] 19.06169

Самая большая ошибка.

Вообще, в данном случае SSA сработал не очень хорошо, так как последний (отрезанный) период повторяет предыдущий, а SSA пытается учесть предыдущую информацию.

Итоги

В отчете мы сравнили три подхода к построению прогноза для временных рядов: ETS, SARIMA, R-SSA, V-SSA. Среднеквадратичные ошибки:

all.mse <- c(ETS = jj.ets.mse, ETS_DUMPED = jj.ets.d.mse, SARIMA = jj.arima.mse, SARIMA_BOX_COX = jj.log.arima.mse, R_SSA = jj.rssa.mse, V_SSA = jj.vssa.mse)
all.mse
##            ETS     ETS_DUMPED         SARIMA SARIMA_BOX_COX          R_SSA 
##      0.4987922      0.1765835      0.2710524      0.1453750      1.2799866 
##          V_SSA 
##     19.0616919

Наименьшей ошибки удалось добиться с помощью SARIMA, когда к исходным данным было применено преобразование Бокса-Кокса. Эта ошибка сравнима для ETS с демпированным трендом.

На наш взгляд, ETS для данного случая подходит лучше потому, что этот метод позволяет учесть мультипликативность модели, тогда как для SARIMA пришлось делать преобразование Бокса-Кокса.

Прогнозы с помощью SSA показали худшие результаты. Возможно, это связано со структурой ряда ближе к концу. Давайте сделаем предсказания для середины ряда.

Дополнение. Предсказание для середины ряда

Мы возьмем первые \(40\) наблюдений для оценки параметров моделей и, как и раньше, сделаем предсказание на год вперед.

mid.split <- train.test.split(head(JohnsonJohnson, 44), test.size)

SSA

Определим количество компонент:

jj.mid.ssa <- ssa(mid.split$Train)
plot(wcor(jj.mid.ssa))

Нужно взять первые четыре компоненты для восстановления сигнала. Предскажем значения:

jj.mid.rssa.forecast <- forecast(jj.mid.ssa, groups = list(Signal = c(1:4)), len = test.size,
                                 method = "recurrent", interval = "prediction", only.intervals = T, level = 0.99)
plot(jj.mid.rssa.forecast)
lines(JohnsonJohnson)

Посмотрим поближе:

tser.my.plot(mid.split$Test, jj.mid.rssa.forecast$mean, component.labs = c('Predicted'), plot.resid = F)

В целом предсказание повторяет форму ряда, хотя и заметно отклоняется от истинных значений.

jj.mid.rssa.mse <- mean((mid.split$Test - jj.mid.rssa.forecast$mean) ^ 2)
jj.mid.rssa.mse
## [1] 0.2899964

ARIMA

jj.mid.arima <- auto.arima(mid.split$Train, d = 1, D = 1, ic = "aicc", stepwise = F, trace = F)
jj.mid.arima
## Series: mid.split$Train 
## ARIMA(1,1,0)(0,1,1)[4] 
## 
## Coefficients:
##           ar1     sma1
##       -0.5095  -0.4913
## s.e.   0.1634   0.1580
## 
## sigma^2 estimated as 0.01845:  log likelihood=20.57
## AIC=-35.14   AICc=-34.36   BIC=-30.47
jj.mid.arima.forecast <- forecast(jj.mid.arima, h = test.size)
plot(jj.mid.arima.forecast)
lines(JohnsonJohnson)

tser.my.plot(mid.split$Test, jj.mid.arima.forecast$mean, component.labs = c('Predicted'), plot.resid = F)

Предсказанные значения повторяют форму, но сильно отклоняются от значений ряда.

jj.mid.arima.mse <- mean((mid.split$Test - jj.mid.arima.forecast$mean) ^ 2)
jj.mid.arima.mse
## [1] 0.7418953

ETS

jj.mid.ets.d <- ets(mid.split$Train, ic = "aicc", opt.crit = "lik", damped = T)
jj.mid.ets.d
## ETS(M,Ad,A) 
## 
## Call:
##  ets(y = mid.split$Train, damped = T, opt.crit = "lik", ic = "aicc") 
## 
##   Smoothing parameters:
##     alpha = 0.1254 
##     beta  = 0.0902 
##     gamma = 0.0312 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 0.6419 
##     b = -0.0069 
##     s = -0.1466 0.1689 -0.0046 -0.0177
## 
##   sigma:  0.1018
## 
##         AIC        AICc         BIC 
## -16.0417371  -8.4555302   0.8470574
jj.mid.ets.d.forecast <- forecast(jj.mid.ets.d, h = test.size)
plot(jj.mid.ets.d.forecast)
lines(JohnsonJohnson)

tser.my.plot(mid.split$Test, jj.mid.ets.d.forecast$mean, component.labs = c('Predcited'), plot.resid = F)

jj.mid.ets.d.mse <- mean((mid.split$Test - jj.mid.ets.d.forecast$mean) ^ 2)
jj.mid.ets.d.mse
## [1] 0.4281315

Сравнение результатов

c(SSA = jj.mid.rssa.mse, SARIMA = jj.mid.arima.mse, ETS = jj.mid.ets.d.mse)
##       SSA    SARIMA       ETS 
## 0.2899964 0.7418953 0.4281315

Теперь результаты методов отличаются не так сильно. Можно предположить, что к концу ряда структура меняется, из-за этого SSA делает столь неточные прогнозы. arima и ets к изменению поведения не так чувствительны: предсказания сильно зависят только от нескольких последних значений ряда, вклад предыдущих сильно ограничен, поэтому важно, чтобы поведение ближайших точек совпадало с точками предсказания. Это требование как раз выполняется для конца ряда.

Источники

  1. https://otexts.com/fpp2/ – обзор моделей ETS и ARIMA
  2. https://faculty.chicagobooth.edu/jeffrey.russell/teaching/finecon/handouts/notes2.pdf – обзор MA(q) модели
  3. https://faculty.chicagobooth.edu/jeffrey.russell/teaching/finecon/handouts/notes2.pdf – выбор параметров ARMA(p, q)
  4. https://www.statisticshowto.datasciencecentral.com/unit-root/ –краткий обзор тестов на стационарность
  5. https://faculty.chicagobooth.edu/ruey.tsay/teaching/uts/lec11-08.pdf – более делтальный разбор Unit-Root tests
  6. “Метод”Гусеница“-SSA: прогноз временных рядов”, Н.Э. Голяндина