Najmniej głupi sposób na prognozowanie krótkich wielowymiarowych szeregów czasowych

Muszę prognozować następujące 4 zmienne dla 29. jednostki czasu. Mam dane historyczne o wartości około 2 lat, gdzie 1, 14 i 27 to ten sam okres (lub pora roku). Na koniec dokonuję dekompozycji w stylu Oaxaca-Blindera na , , i p .w d w c $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

Uważam, że $W$ można oszacować za pomocą $p\cdot wd + (1 - p)\cdot wc$ plus błąd pomiaru, ale widać, że $W$ zawsze znacznie przekracza tę ilość z powodu marnotrawstwa, błędu przybliżenia lub kradzieży.

Oto moje 2 pytania.

1. Moją pierwszą myślą było wypróbowanie wektorowej autoregresji tych zmiennych z 1 opóźnieniem i egzogenną zmienną czasu i okresu, ale wydaje się to złym pomysłem, biorąc pod uwagę, jak mało danych mam. Czy istnieją metody szeregów czasowych, które (1) działają lepiej w obliczu „mikroliczności” i (2) mogłyby wykorzystać związek między zmiennymi?

2. Z drugiej strony, moduły wartości własnych dla VAR są mniejsze niż 1, więc nie sądzę, żebym musiał się martwić o niestacjonarność (chociaż test Dickeya-Fullera sugeruje inaczej). Prognozy wydają się w większości zgodne z prognozami elastycznego modelu jednowymiarowego z trendem czasowym, z wyjątkiem i , które są niższe. Współczynniki opóźnień wydają się w większości rozsądne, choć w większości są nieznaczne. Współczynnik trendu liniowego jest znaczący, podobnie jak niektóre z manekinów okresu. Czy istnieją jednak teoretyczne powody, aby preferować to prostsze podejście niż model VAR?$W$$p$

Pełne ujawnienie: zadałem podobne pytanie w sprawie Statalist bez odpowiedzi.

Rozumiem, że to pytanie spotyka się tutaj od lat, ale nadal przydatne mogą być następujące pomysły:

1. Jeśli istnieją powiązania między zmiennymi (a wzór teoretyczny nie działa tak dobrze), PCA można wykorzystać do systematycznego wyszukiwania zależności (liniowych). Pokażę, że działa to dobrze dla danych w tym pytaniu.

2. Biorąc pod uwagę, że nie ma zbyt wielu danych (łącznie 112 liczb), można oszacować tylko kilka parametrów modelu ( np. Dopasowanie pełnych efektów sezonowych nie jest opcją), a wypróbowanie modelu niestandardowego może mieć sens.

Oto jak stworzyłbym prognozę zgodnie z następującymi zasadami:

Krok 1. Możemy użyć PCA do ujawnienia zależności w danych. Używając R, z danymi przechowywanymi w 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

To pokazuje, że dwa pierwsze główne składniki wyjaśniają 97% wariancji, a użycie trzech komponentów pokrywa 99,4% wariancji. Wystarczy więc stworzyć model na pierwsze dwa lub trzy komputery. (Dane w przybliżeniu spełniają .)$W = 0.234\, wd - 1.152\, wc - 8.842 \,p$

Wykonanie PCA wymagało znalezienia macierzy ortogonalnej . Przestrzeń takich matryc jest 6-wymiarowa, więc oszacowaliśmy 6 parametrów. (Ponieważ tak naprawdę używamy tylko PC1 poniżej, może to być mniej „efektywnych” parametrów).$4\times 4$

Krok 2. W PC1 jest wyraźny trend:

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

Tworzę kopię wyników na PC po usunięciu tego trendu:

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

Rysowanie wyników innych komputerów nie ujawnia żadnych wyraźnych trendów, więc pozostawiam je bez zmian.

Ponieważ wyniki PC są wyśrodkowane, trend przechodzi przez środek masy próbki PC1, a dopasowanie trendu odpowiada tylko oszacowaniu jednego parametru.

Krok 3. Para wykresu rozrzutu nie wykazuje wyraźnej struktury, więc modeluję komputery jako niezależne:

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

Krok 4. W PC1 występuje wyraźna okresowość, z opóźnieniem 13 (jak sugeruje pytanie). Można to postrzegać na różne sposoby. Na przykład autokorelacja opóźnienia 13 wykazuje się jako istotnie różna od 0 w korelogramie:

> acf(y2[,1])

(Okresowość jest wizualnie bardziej uderzająca przy drukowaniu danych wraz z przesuniętą kopią).

Ponieważ chcemy utrzymać liczbę oszacowanych parametrów na niskim poziomie, a ponieważ korelogram pokazuje opóźnienie 13 jako jedyne opóźnienie ze znaczącym udziałem, modeluję PC1 jako , gdzie są niezależne i standardowo rozkładają się normalnie (tj. Jest to proces AR (13) z większością współczynników ustalonych na 0). Łatwym sposobem oszacowania i jest użycie funkcji:$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

Jako test wiarygodności wykreślam dane (czarne) wraz z losową trajektorią mojego modelu dla PC1 (niebieski), począwszy od roku:

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

Niebieski, symulowany fragment ścieżki wygląda jak rozsądna kontynuacja danych. Korelogramy dla PC2 i PC3 nie wykazują znaczących korelacji, więc modeluję te elementy jako biały szum. PC4 wykazuje korelacje, ale tak mało przyczynia się do całkowitej wariancji, że wydaje się, że nie warto modelować, a ja również modeluję ten komponent jako biały szum.

Tutaj umieściliśmy dwa dodatkowe parametry. Daje nam to w sumie dziewięć parametrów w modelu (w tym PCA), co nie wydaje się absurdalne, gdy zaczynamy od danych składających się ze 112 liczb.

Prognoza. Możemy uzyskać prognozę numeryczną, pomijając szum (aby uzyskać średnią) i odwracając 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

Pasma niepewności można uzyskać analitycznie lub po prostu za pomocą Monte Carlo:

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$