Наименее глупый способ прогнозирования коротких многомерных временных рядов

Мне нужно спрогнозировать следующие 4 переменные для 29-й единицы времени. У меня есть исторические данные примерно за 2 года, где 1 и 14 и 27 - все один и тот же период (или время года). В конце я делаю разложение в стиле Оахака-Блиндера на $W$ , $wd$ , $wc$ и $p$ .

time    W               wd              wc               p
1       4.920725        4.684342        4.065288        .5962985
2       4.956172        4.73998         4.092179        .6151785
3       4.85532         4.725982        4.002519        .6028712
4       4.754887        4.674568        3.988028        .5943888
5       4.862039        4.758899        4.045568        .5925704
6       5.039032        4.791101        4.071131        .590314
7       4.612594        4.656253        4.136271        .529247
8       4.722339        4.631588        3.994956        .5801989
9       4.679251        4.647347        3.954906        .5832723
10      4.736177        4.679152        3.974465        .5843731
11      4.738954        4.759482        4.037036        .5868722
12      4.571325        4.707446        4.110281        .556147
13      4.883891        4.750031        4.168203        .602057
14      4.652408        4.703114        4.042872        .6059471
15      4.677363        4.744875        4.232081        .5672519
16      4.695732        4.614248        3.998735        .5838578
17      4.633575        4.6025          3.943488        .5914644
18      4.61025         4.67733         4.066427        .548952
19      4.678374        4.741046        4.060458        .5416393
20      4.48309         4.609238        4.000201        .5372143
21      4.477549        4.583907        3.94821         .5515663
22      4.555191        4.627404        3.93675         .5542806
23      4.508585        4.595927        3.881685        .5572687
24      4.467037        4.619762        3.909551        .5645944
25      4.326283        4.544351        3.877583        .5738906
26      4.672741        4.599463        3.953772        .5769604
27      4.53551         4.506167        3.808779        .5831352
28      4.528004        4.622972        3.90481         .5968299

Я считаю, что $W$ можно аппроксимировать как плюс ошибка измерения, но вы можете видеть, что всегда значительно превышает это количество из-за потерь, ошибки аппроксимации или кражи.$p\cdot wd + (1 - p)\cdot wc$$W$

Вот мои 2 вопроса.

1. Моей первой мыслью было попробовать векторную авторегрессию для этих переменных с 1 лагом и экзогенной переменной времени и периода, но это кажется плохой идеей, учитывая, как мало у меня данных. Существуют ли какие-либо методы временных рядов, которые (1) работают лучше перед лицом "микросчетности" и (2) могли бы использовать связь между переменными?

2. С другой стороны, все модули собственных значений для VAR меньше 1, поэтому я не думаю, что мне нужно беспокоиться о нестационарности (хотя тест Дики-Фуллера предполагает иное). Прогнозы в основном соответствуют прогнозам гибкой однофакторной модели с временным трендом, за исключением и , которые ниже. Коэффициенты на лагах кажутся в основном разумными, хотя по большей части они незначительны. Коэффициент линейного тренда является значительным, как и некоторые манекены периода. Тем не менее, есть ли теоретические причины, чтобы предпочесть этот более простой подход модели VAR?$W$$p$

Полное раскрытие: я задал аналогичный вопрос на Statalist без ответа.

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

1. Если между переменными есть связи (а теоретическая формула работает не очень хорошо), PCA можно использовать для систематического поиска (линейных) зависимостей. Я покажу, что это хорошо работает для данных в этом вопросе.

2. Учитывая, что данных не так много (всего 112 номеров), можно оценить только несколько параметров модели ( например, подбор полных сезонных эффектов не представляется возможным), и может иметь смысл попробовать пользовательскую модель.

Вот как я бы сделал прогноз, следуя этим принципам:

Шаг 1. Мы можем использовать PCA для выявления зависимостей в данных. Используя R, данные хранятся в x:

> library(jvcoords)
> m <- PCA(x)
> m
PCA: mapping p = 4 coordinates to q = 4 coordinates

                              PC1         PC2          PC3          PC4
standard deviation     0.18609759 0.079351671 0.0305622047 0.0155353709
variance               0.03463231 0.006296688 0.0009340484 0.0002413477
cum. variance fraction 0.82253436 0.972083769 0.9942678731 1.0000000000

$W = 0.234\, wd - 1.152\, wc - 8.842 \,p$

$4\times 4$

Шаг 2. Есть четкая тенденция в ПК1:

> t <- 1:28
> plot(m$y[,1], type = "b", ylab = "PC1")
> trend <- lm(m$y[,1] ~ t)
> abline(trend)

Я создаю копию результатов ПК, убрав эту тенденцию:

> y2 <- m$y
> y2[,1] <- y2[,1] - fitted(trend)

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

Поскольку оценки ПК центрированы, тренд проходит через центр масс образца ПК1, и подгонка тренда соответствует оценке только одного параметра.

Шаг 3. График разброса пар не показывает четкой структуры, поэтому я моделирую ПК как независимые:

> pairs(y2, asp = 1, oma = c(1.7, 1.7, 1.7, 1.7))

Шаг 4. Существует четкая периодичность в ПК1 с лагом 13 (как предполагает вопрос). Это можно увидеть по-разному. Например, автокорреляция с лагом 13 значительно отличается от 0 на коррелограмме:

> acf(y2[,1])

(Периодичность визуально более поразительна при нанесении данных вместе со смещенной копией.)

$y^{(1)}_{t+13} = \alpha_{13} y^{(1)}_t + \sigma \varepsilon_{t+13}$$\varepsilon_t$$\alpha_{13}$$\sigma$lm()

> lag13 <- lm(y2[14:28,1] ~ y2[1:15,1] + 0)
> lag13

Call:
lm(formula = y2[14:28, 1] ~ y2[1:15, 1] + 0)

Coefficients:
y2[1:15, 1]  
     0.6479  

> a13 <- coef(lag13)
> s13 <- summary(lag13)$sigma

В качестве теста правдоподобия я нанесу данные (черный) вместе со случайной траекторией моей модели для ПК1 (синий), рассчитанной на один год в будущем:

t.f <- 29:41
pc1 <- m$y[,1]
pc1.f <- (predict(trend, newdata = data.frame(t = t.f))
          + a13 * y2[16:28, 1]
          + rnorm(13, sd = s13))
plot(t, pc1, xlim = range(t, t.f), ylim = range(pc1, pc1.f),
     type = "b", ylab = "PC1")
points(t.f, pc1.f, col = "blue", type = "b")

Синий, смоделированный участок пути выглядит как разумное продолжение данных. Коррелограммы для ПК2 и ПК3 не показывают значительных корреляций, поэтому я моделирую эти компоненты как белый шум. PC4 действительно показывает корреляции, но вносит настолько малый вклад в общую дисперсию, что кажется не стоит моделировать, и я также моделирую этот компонент как белый шум.

Здесь мы установили еще два параметра. Это подводит нас к девяти параметрам в модели (включая PCA), что не кажется абсурдным, когда мы начали с данных, состоящих из 112 чисел.

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

> pc1.f <- predict(trend, newdata = data.frame(t = t.f)) + a13 * y2[16:28, 1]
> y.f <- data.frame(PC1 = pc1.f, PC2 = 0, PC3 = 0, PC4 = 0)
> x.f <- fromCoords(m, y.f)
> rownames(x.f) <- t.f
> x.f
          W       wd       wc         p
29 4.456825 4.582231 3.919151 0.5616497
30 4.407551 4.563510 3.899012 0.5582053
31 4.427701 4.571166 3.907248 0.5596139
32 4.466062 4.585740 3.922927 0.5622955
33 4.327391 4.533055 3.866250 0.5526018
34 4.304330 4.524294 3.856824 0.5509898
35 4.342835 4.538923 3.872562 0.5536814
36 4.297404 4.521663 3.853993 0.5505056
37 4.281638 4.515673 3.847549 0.5494035
38 4.186515 4.479533 3.808671 0.5427540
39 4.377147 4.551959 3.886586 0.5560799
40 4.257569 4.506528 3.837712 0.5477210
41 4.289875 4.518802 3.850916 0.5499793

Полосы неопределенности могут быть получены аналитически или просто с использованием метода Монте-Карло:

N <- 1000 # number of Monte Carlo samples
W.f <- matrix(NA, N, 13)
for (i in 1:N) {
    y.f <- data.frame(PC1 = (predict(trend, newdata = data.frame(t = t.f))
              + a13 * y2[16:28, 1]
              + rnorm(13, sd = s13)),
              PC2 = rnorm(13, sd = sd(y2[,2])),
              PC3 = rnorm(13, sd = sd(y2[, 3])),
              PC4 = rnorm(13, sd = sd(y2[, 4])))
    x.f <- fromCoords(m, y.f)
    W.f[i,] <- x.f[, 1]
}
bands <- apply(W.f, 2,
               function(x) quantile(x, c(0.025, 0.15, 0.5, 0.85, 0.975)))
plot(t, x$W, xlim = range(t, t.f), ylim = range(x$W, bands),
     type = "b", ylab = "W")
for (b in 1:5) {
    lines(c(28, t.f), c(x$W[28], bands[b,]), col = "grey")
}

$W$