Как найти изменяющиеся во времени коэффициенты для модели VAR с помощью фильтра Калмана - PullRequest
Новая
       58

Как найти изменяющиеся во времени коэффициенты для модели VAR с помощью фильтра Калмана

1 голос
/ 06 апреля 2019

Я пытаюсь написать код на R, чтобы воспроизвести модель, которую я нашел в этой статье .

Идея состоит в том, чтобы смоделировать сигнал как модель VAR, но подгонятькоэффициенты по модели фильтра Калмана.Это, по сути, позволило бы мне создать надежную изменяющуюся во времени модель VAR (p) и до некоторой степени анализировать нестационарные данные.

Модель для отслеживания коэффициентов:

X (t) = F (t) X (t− 1) + W (t)

Y (t)= H (t) X (t) + E (t),

, где H (t) - произведение Кронекера между запаздывающими измерениями в моем временном ряду Y и единичным вектором,и X (t) выполняет роль коэффициентов регрессии. F (t) принимается за единичную матрицу, поскольку это должно означать, что мы предполагаем, что коэффициенты будут развиваться как случайное блуждание.

В статье из W (T) ковариационная матрица шума состояний Q (t) сначала выбирается при 10 ^ -3, а затем устанавливается на основепо какой-то итерационной схеме.Начиная с E (t) ковариационная матрица шума состояний равна R (t) , заменяемой ковариацией шумового члена, необъясненного моделью: Y (t) - H (t) Xhat (t)

У меня есть априорная ковариационная матрица ошибки оценки (обозначается в статье Σ ), записанная как P (на основе других источников) и a posteriori как Pmin , так как он будет использоваться в следующей рекурсии как a priori , если это делаетсмысл.

Пока что я написал следующее, основываясь на статьях Приложение A 1.2 

Y <-                                  *my timeseries, for test purposes two channels of 3000 points*
F <- diag(8)                          # F is (m^2*p by m^2 *p) where m=2 dimensions and p =2 lags
H <- diag(2) %x% t(vec(Y[,1:2]))      #the Kronecker product of vectorized lags Y-1 and Y-2
Xhatminus <- matrix(1,8,1)            # an arbitrary *a priori* coefficient matrix
Q <- diag(8)%x%(10**-7)               #just a number a really low number diagonal matrix, found it used in some examples
R<- 1                                 #Didnt know what else to put here just yet
Pmin = diag(8)                        #*a priori* error estimate, just some 1-s...

Теперь следует начать рекурсию.Для тестирования я просто взял первые 3000 баллов за один тест моих данных.

   Xhatstorage <- matrix(0,8,3000)

   for(j in 3:3000){
            H <- diag(2) %x% t(vec(Y[,(j-2):(j-1)]))
            K <- (Pmin %*% t(H)) %*% solve((H%*%Pmin%*%t(H) + R)) ##Solve gives inverse matrix ()^-1
            P <- Pmin - K%*% H %*% Pmin
            Xhatplus <- F%*%( Xhatminus + K%*%(Y[,j]-H%*%Xhatminus) )
            Pplus <- (F%*% P %*% F)  + Q
            Xhatminus <- Xhatplus
            Xhatstorage[,j] <- Xhatplus 
            Pmin <- Pplus
    }

Я извлек Xhatplus значения в матрицу хранения и использовал их для записи с ними этой примитивной модели VAR: 

Yhat<-array(0,3000)
for(t in 3:3000){
Yhat[t]<- t(vec(Y[,(t-2)])) %*% Xhatstorage[c(1,3),t] + t(vec(Y[,(t-1)])) %*% Xhatstorage[c(2,4),t]
}

Похоже на это .

Синяя линия - это VAR с коэффициентами, найденными фильтром Калмана, черный - исходные данные ..

У меня возникли проблемы с пониманием того, как я могу лучше оценить свои коэффициенты?Почему это так?

Как мне лучше выбрать первые априорные и апостериорные оценки, чтобы начать рекурсию?В настоящее время, я уверен, что добавление большего количества лагов в VAR не проблема, я не знаю, как выбрать начальные значения для Pmin и Xhatmin.В большинстве мест я собирал это с самого начала из произвольных 0 предположений в игрушечных моделях, но в этом случае выбор любой из упомянутых матриц в качестве 0 просто свернет весь алгоритм.

Наконец, является ли эта рекурсия даже правильной реализацией Oya et al , описанной в статье?Я знаю, что мне все еще не хватает оценки R, основанной на ранее необъяснимых ошибках (V (t) в Приложении A 1.2), но в целом?

...