Данные

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

library(Rssa)
library(pracma)
library(zoo)
library(stats)
library(plotrix)
library(MASS)

#setwd("/media/ekaterina/Data/Study/Magistracy/10 term/Time series (pr)/SSA_real/")
setwd("E:/Study/Magistracy/10 term/Time series (pr)/Reports")
data_raw <- read.csv("milk.csv")
data_zoo <- zoo(data_raw[,2], as.Date(data_raw$Month, format = '%Y/%m/%d'))
data <- ts(data_zoo)
plot(data_zoo)

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

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

По пикам периодограммы можно говорить о наличии сезонных компонент с частотами 1/12, 2/12, 3/12, 4/12, 5/12 и 6/12.

SSA

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

Выделение тренда

Будем считать, что тренд сложной формы. Для выделения тренда выберем длину небольшую длину окна, кратную периоду - 24.

ssa1 <- ssa(data, L = 24)
plot(ssa1)

plot(ssa1, type = "vectors", idx = 1:15)

plot(ssa1, type = "paired", idx = 1:15)

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

plot(ssa1, type = "series", groups = as.list(1:15))

plot(wcor(ssa1))

К тренду отнесем медленно меняющиеся компоненты. 1 и 8.

recon1 <- reconstruct(ssa1, groups = list(c(1,8)))
trend <- recon1$F1
plot(data, type = "l")
lines(trend, col = "red")

Сезонность

resid1 <- residuals(recon1)
plot(resid1, type= "l", main = "Остатки")

Для выделения сезонности выберем окно кратное периоду и близкое к половине ряда.

ssa2 <- ssa(resid1)

plot(ssa2)

plot(ssa2, type = "vectors", idx = 1:15)

plot(ssa2, type = "paired", idx = 1:15)

plot(ssa2, type = "series", groups = as.list(1:15))

#plot(wcor(ssa2)) # w-correlation matrix plot

На графиках пар собственных векторов видно, что компоненты, которые описывают колебания с частотами 4/12 и 3/12, смешались и выделились плохо. Это также видно на матрице взешенных корреляций.

plot(wcor(ssa2))

На этой матрице видно, что 1-2, 3-4, 5-6 и 11 компоненты отделились хорошо, однако 7-8 и 9-10 смешались.

Теперь оценим частоты синусов, которые описавают с 1 по 11 компоненты.

parestimate(ssa2, list(1:11), method =  "esprit-ls")
##    period     rate   |    Mod     Arg  |     Re        Im
##     3.002   0.002978 |  1.00298   2.09 | -0.50010   0.86941
##    -3.002   0.002978 |  1.00298  -2.09 | -0.50010  -0.86941
##     2.400   0.002652 |  1.00266   2.62 | -0.86851   0.50100
##    -2.400   0.002652 |  1.00266  -2.62 | -0.86851  -0.50100
##     2.000   0.001970 |  1.00197   3.14 | -1.00197   0.00000
##     3.987   0.001552 |  1.00155   1.58 | -0.00508   1.00154
##    -3.987   0.001552 |  1.00155  -1.58 | -0.00508  -1.00154
##    12.032   0.000211 |  1.00021   0.52 |  0.86690   0.49890
##   -12.032   0.000211 |  1.00021  -0.52 |  0.86690  -0.49890
##     6.002  -0.001528 |  0.99847   1.05 |  0.49951   0.86454
##    -6.002  -0.001528 |  0.99847  -1.05 |  0.49951  -0.86454

То есть, 1-11 компоенты описыют все частоты, которые мы выделили на периодограмме. Будем их использовать для выделения сезонности.

recon2 <- reconstruct(ssa2, groups=list(c(1:11)))
seasonality <- recon2$F1
resid2 <- residuals(recon2)
plot(recon2)

Остатки

Посмотрим на график остатков.

plot(resid2, main = "residuals")

acf(resid2)

По acf остатков можно предположить, что это белый шум.

Filter-adjusted O-SSA

Так как есть компоненты, которые плохо отделились, попробуем их отделить получше. Будем считать, что у нас есть сигнал без шума.

signal <- trend + seasonality
plot(signal, type = "l")

signal.ssa <- ssa(signal)

plot(signal.ssa)

plot(signal.ssa, type = "vectors", idx = 1:15)

plot(signal.ssa, type = "paired", idx = 1:15)

plot(signal.ssa, type = "series", groups = as.list(1:15))

plot(wcor(signal.ssa, groups = 1:14))

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

Видим, что смешались компоненты тренда 8 и 9, И компоненты сезонности 10 и 11, 12 и 13.

fs <- fossa(signal.ssa, nested.groups = list(8:9, 10:11, 12:13))

plot(fs)

plot(fs, type = "vectors", idx = 1:15)

plot(fs, type = "paired", idx = 1:15)

plot(fs, type = "series", groups = as.list(1:15))

plot(wcor(fs, groups = 1:14))

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

recon3 <- reconstruct(fs, groups=list(c(1,12,13), c(2:11, 14)))
plot(recon3)

trend.fs <- recon3$F1
seasonalty.fs <- recon3$F2

Корни характерестического полинома

Сигнальные корни, это те корни, у которых модуль для прямого или развернутого ряда больше 1.

signal <- trend + seasonality
plot(signal, type = "l")

ssa_signal <- ssa(signal)
n_eig <- list(1:14)

l1<- lrr(ssa_signal, group = n_eig)
r1 <- Rssa::roots(l1)
main.roots <- r1[Mod(r1) > 1]

l2 <- lrr(ssa_signal, reverse = TRUE, group = n_eig)
r2 <- Rssa::roots(l2)
main.roots <- c(main.roots, 1/r2[Mod(r2) > 1])
main.roots
##  [1]  1.0015702+0.0914841i  1.0015702-0.0914841i -1.0038726+0.0000000i
##  [4] -0.4995689+0.8703577i -0.4995689-0.8703577i -0.8687201+0.5010459i
##  [7] -0.8687201-0.5010459i -0.0027539+1.0022842i -0.0027539-1.0022842i
## [10]  1.0018428+0.0000000i  0.8671275+0.4990176i  0.8671275-0.4990176i
## [13]  0.4989110-0.8643013i  0.4989110+0.8643013i
plot(main.roots, col = "blue")
circle(0, 0, 1, add = TRUE)

2*pi/Arg(main.roots)
##  [1]  68.979348 -68.979348   2.000000   3.003629  -3.003629   2.399596
##  [7]  -2.399596   3.993016  -3.993016        Inf  12.032198 -12.032198
## [13]  -5.999535   5.999535

Видим, что оценки частот для главных корней совпадают с нашими ожиданиями. 11 корней для сезонности и 3 корня для тренда.

Cadzow

Попробуем с помщью метода Cadzow получить ряд конечного ранга для наших данных. Выберем ранг равным 13. 11 компонент для описания сезонности и 2 компоненты для описания тренда.

ssa.dat <- ssa(data)
cad <- cadzow(ssa.dat, rank = 13, correct = FALSE)
library(lattice)
xyplot(cbind(Original = data, Cadzow = cad), superpose = TRUE)

ssa.cad <- ssa(cad)
plot(signal.ssa)

Видим, что после 13 собственные числа начинают резко убывать к нулю.

Сравним Cadzow и SSA восстановление

rec.ssa <- reconstruct(ssa.dat, groups = list(1:13))
rec.cad <- reconstruct(ssa.cad, groups = list(1:13))

xyplot(cbind(Original = data, Cadzow = rec.cad$F1, SSA = rec.ssa$F1), superpose = TRUE)

В целом результаты выглядят очень схожими.

Прогноз

В SSA используется два типа прогноза - реккурентный и векторный.

Разделим выборку на тестовую и тренировачную. Будем предсказывать на 3 периода.

Сравним 2 периода предсказания с истинными данными и оценим MSE для каждого из прогнозов.

Реккурентный

h = 36
train <- data[1:(length(data)-h)]
ssa.train <- ssa(train)

pred.ssa.rec <- predict(ssa.train, group = list(1:14), method = "bootstrap-recurrent", len = h)
start.point <- length(data) - 24

plot(data, type = "l", xlim = c(1,(start.point + h)), ylim = c(550, 1000))
lines(start.point:(start.point + h-1),pred.ssa.rec[,1], type = "l", col = "red")
lines(start.point:(start.point + h-1),pred.ssa.rec[,2], type = "l", col = "blue", lty = 2)
lines(start.point:(start.point + h-1),pred.ssa.rec[,3], type = "l", col = "blue", lty = 2)

mse.rec <- mean((data[(length(data) - h+1): length(data)] - pred.ssa.rec[,1])^2)

Для реккурентного прогноза MSE = 1646.82.

Векторный

pred.ssa.vec <- predict(ssa.train, group = list(1:14), method = "bootstrap-vector", len = h)

start.point <- length(data) - 24

plot(data, type = "l", xlim = c(1,(start.point + h)), ylim = c(550, 1000))
lines(start.point:(start.point + h-1),pred.ssa.vec[,1], type = "l", col = "red")
lines(start.point:(start.point + h-1),pred.ssa.vec[,2], type = "l", col = "blue", lty = 2)
lines(start.point:(start.point + h-1),pred.ssa.vec[,3], type = "l", col = "blue", lty = 2)

mse.vec <- mean((data[(length(data) - h+1): length(data)] - pred.ssa.vec[,1])^2)

Для векторного прогноза MSE = 1810.18. Для реккурентного прогноза ошибка меньше.

Автоматическая идентификация

Частотная группировка

gr <- grouping.auto(ssa1, grouping.method = "pgram", groups = c(1,8), base = "series", treshold = 0.8)
plot(reconstruct(ssa1, groups = gr))

Кластерная группировка

gr <- grouping.auto(ssa1, grouping.method = "wcor", groups =c(1,8) ,nclust = 1)
plot(reconstruct(ssa1, groups = gr))

Заполнение пропусков

Subspace based

Добавим пропуски к данным

data.na <- data
data.na[50:60] <- NA
plot(data.na)

na.ssa <- ssa(data.na)
## Warning: Some field elements were not covered by shaped window. 49 elements
## will be ommited
g <- gapfill(na.ssa, groups = list(1:14))
plot(g)
lines(data, col = "blue")

mse.g <- mean((g[50:60] - data[50:60])^2)
mse.g
## [1] 120.3642

Интервальный

g <- igapfill(na.ssa, groups = list(1:14))
plot(g)
lines(data, col = "blue")

mse.ig <- mean((g[50:60] - data[50:60])^2)
mse.ig
## [1] 114.5912

Ошибка немного меньше для интервального заполнения.

Projecttion SSA

Projecttion SSA применяется для улучшения разделимости линейного тренда и сезонности.

Смоделиреум данные с линейным трендом.

x <- 1:100
data.l <- x*2 + 3 + 100*cos(2*pi*(1/25)*x) + cos(2*pi*(1/5)*x) + rnorm(50)
plot(data.l, type = "l")

Применим SSA с двойным центрированием для лучшего отделения линейного тренда.

ssa.l <- ssa(data.l, L = 25, column.projector='centering', row.projector='centering')


plot(ssa.l, type = "vectors") # Eigenvectors

plot(ssa.l, type = "paired") # Pairs of eigenvectors

plot(ssa.l, type = "series", groups = as.list(1:10))

rec.l <- reconstruct(ssa.l, groups = list(1:2))
trend.l <- rec.l$F1

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

Теплицев SSA

Теплицев SSA применяется как улучшение basic-SSA для стационарных рядов.

Смоделируем стационарный ряд

x <- 1:100
data.s <- 5*cos(2*pi*(1/25)*x) + cos(2*pi*(1/5)*x) + rnorm(2)
plot(data.s, type = "l")

ssa.s <- ssa(data.s, kind = "toeplitz-ssa")

plot(ssa.s, type = "paired") # Pairs of eigenvectors

rec.s <- reconstruct(ssa.s, groups = list(c(1,2,4,5)))
trend.s <- rec.s$F1
resid.s <- residuals(rec.s)

plot(data.s, type = "l")
lines(trend.s, col = "red")