Mam model mieszanego efektu (w rzeczywistości uogólniony model mieszany dodatku), który daje mi prognozy dla szeregów czasowych. Aby przeciwdziałać autokorelacji, używam modelu corCAR1, biorąc pod uwagę fakt, że brakuje mi danych. Dane powinny dać mi całkowite obciążenie, więc muszę sumować przez cały przedział prognozowania. Ale powinienem również uzyskać oszacowanie błędu standardowego tego całkowitego obciążenia.
Jeśli wszystkie prognozy byłyby niezależne, można to łatwo rozwiązać przez:
z V a r ( E [ X i ] ) = S E ( E [ X i ] ) 2
Problem polega na tym, że przewidywane wartości pochodzą z modelu, a oryginalne dane mają autokorelację. Cały problem prowadzi do następujących pytań:
- Czy mam rację, zakładając, że SE obliczonych prognoz może być interpretowane jako źródło wariancji oczekiwanej wartości tej prognozy? Mam tendencję do interpretowania prognoz jako „średnich prognoz”, a zatem sumowania całego zestawu środków.
- Jak włączyć autokorelację do tego problemu lub czy mogę bezpiecznie założyć, że nie wpłynie to zbytnio na wyniki?
Jest to przykład w R. Mój prawdziwy zestaw danych zawiera około 34 000 pomiarów, więc problemem jest skalowalność. To jest powód, dla którego modeluję autokorelację w ciągu każdego miesiąca, w przeciwnym razie obliczenia nie będą już możliwe. To nie jest najbardziej poprawne rozwiązanie, ale najbardziej poprawne nie jest możliwe.
set.seed(12)
require(mgcv)
Data <- data.frame(
dates = seq(as.Date("2011-1-1"),as.Date("2011-12-31"),by="day")
)
Data <- within(Data,{
X <- abs(rnorm(nrow(Data),3))
Y <- 2*X + X^2 + scale(Data$dates)^2
month <- as.POSIXlt(dates)$mon+1
mday <- as.POSIXlt(dates)$mday
})
model <- gamm(Y~s(X)+s(as.numeric(dates)),correlation=corCAR1(form=~mday|month),data=Data)
preds <- predict(model$gam,se=T)
Total <- sum(preds$fit)
Edytować :
Lekcja do nauczenia: najpierw przejrzyj wszystkie próbki we wszystkich plikach pomocy, zanim zaczniesz panikować. W plikach pomocy prognozy.gam mogę znaleźć:
#########################################################
## now get variance of sum of predictions using lpmatrix
#########################################################
Xp <- predict(b,newd,type="lpmatrix")
## Xp %*% coef(b) yields vector of predictions
a <- rep(1,31)
Xs <- t(a) %*% Xp ## Xs %*% coef(b) gives sum of predictions
var.sum <- Xs %*% b$Vp %*% t(Xs)
Co wydaje się być bliskie temu, co chcę zrobić. To wciąż nie mówi mi dokładnie, jak to się robi. Mógłbym dojść do tego, że opiera się on na liniowej macierzy predyktorów. Wszelkie spostrzeżenia są nadal mile widziane.
źródło
Odpowiedzi:
W notacji macierzowej model mieszany można przedstawić jako
y = X * beta + Z * u + epsilon
gdzie X i Z są znanymi macierzami projektowymi dotyczącymi odpowiednio ustalonych efektów i obserwacji efektów losowych.
Zastosowałbym prostą i adekwatną (ale nie najlepszą) transformację do korekcji autokorelacji, która obejmuje utratę pierwszej obserwacji i zastąpienie wektora kolumny [y1, y2, ... yn] mniejszym przez jeden wektor kolumny obserwacyjnej, a mianowicie: [y2 - rho * y1, y3 - rho * y2, ..., yn - rho * y (n-1)], gdzie rho jest oszacowaną wartością seryjnej autokorelacji.
Można tego dokonać mnożąc przez macierz T, tworząc T * y, gdzie pierwszy rząd T składa się w następujący sposób: [-rho, 1, 0, 0, ....], drugi rząd: [0, -rho, 1, 0, 0, ...] itd. Podobnie, inne macierze projektowe są zmieniane na T * X i T * Z. Również zmieniona jest macierz wariancji-kowariancji składników błędów, teraz z niezależnymi warunkami błędu.
Teraz wystarczy obliczyć rozwiązanie za pomocą nowych matryc projektowych.
źródło