Модель VAR в пространстве состояний с фильтрами Калмана - оценка параметров модели и значения

Question

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

В настоящее время у меня есть следующая формулировка двумерной модели p-oder VAR в пространстве состояний:

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

По сути, это представление двумерной модели VAR в пространстве состояний (для двух временных рядов Y1 и Y2):

с коэффициентами A, все найденными в векторе X. Я оцениваю эту модель VAR с помощью следующих итераций фильтра Калмана:

Чтобы запустить алгоритм, мне нужно определить шум процесса Q(t), ковариацию ошибки измерения R(t), априорную ковариационную матрицу ошибки оценки P_{1|0} и начальную оценку X(0).

Я выбираю начальные оценки коэффициентов равными 1, что означает, что X(1) представляет собой (m^2p на 1) матрицу единиц, которая сходится к правильным оценкам в ходе итераций.

Однако определение других начальных параметров сбивает с толку. Так как я хочу получить коэффициенты для модели VAR, "измерения" в моем контексте полностью точны, но я не уверен, означает ли это, что мне нужно взять E(t)=Y(t)-H(t)X(t) или E=0.

В обычном применении фильтра Калмана для целей отслеживания E (t) и через него R (t) представляют ковариацию измерительного устройства. В моем случае я хочу имитировать временной ряд и значения для входных данных " измерения "точно такие, как я хочу, без дальнейших ошибок. Я понимаю, что это означает, что я должен установить E(t)=Y(t)-H(t)X(t) и R(t) = E(t)E(t)^T, и он обновит каждую итерацию.

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

Шум процесса мне более понятен, так как в обычных приложениях он задается аналогично, но не менее сложен в реализации. Я хочу выбрать Q(t), который поддерживает согласованность моей модели, для чего в литературе предлагается использовать статистический тест для нововведений алгоритма фильтра Калмана. Распространенным предложением является индекс нормализованных квадратов инноваций (NIS) (который должен иметь распределение), однако я не смог понять, как реализовать его в моем алгоритме. До сих пор выбор подходил для работы с смоделированными данными (коэффициенты показывают тенденции, которые я встроил в данные).

Наконец, выбирая исходную матрицу ковариации ошибки априорной оценки P. Я начинаю с единичной матрицы размера (m^2p на m^2p), поскольку у меня нет лучшего предположения или понимания, чтобы исправить это.

Мои вопросы, чтобы подвести итог:

1. как лучше всего выбрать член R (t) или обновить E (t) в моем алгоритме?

2. как выбрать Q (t) и, возможно, как реализовать для него тест NIS?

3. на последок как выбрать априори подходящий P?

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

library(fUnitRoots)
library(urca)
library(vars)
library(aod)
library(zoo)
library(tseries)
library(aTSA)

dumdata01<-list(rep(0,times=1000),nrow=(1))
dumdata01<-dumdata01[[1]]

for (i in 2:1000){
  dumdata01[i]<-(dumdata01[i-1]+rnorm(1,mean=0,sd=1))
}
dumdata02<-list(rep(0,times=1000),nrow=(1))
dumdata02<-dumdata02[[1]]
for (i in 2:300){
  dumdata02[i]<-(dumdata02[i-1]+rnorm(1,mean=0,sd=0.5))
  }
for (i in 300:700){
  dumdata02[i]<-(-1)*(dumdata01[i]+(dumdata01[i-1]))/2
}
for (i in 700:850){
  dumdata02[i]<-(dumdata01[i-2]+rnorm(1,mean=0,sd=0.15))
}
for (i in 850:1000){
  dumdata02[i]<-(dumdata02[i-1]+rnorm(1,mean=0,sd=0.5))
}

m<-2
sL <- 1000
Y<-matrix(0,m,sL)
Y[1,]<-dumdata01
Y[2,]<-dumdata02

m<-2
p<-2
F <- diag((m^2)*p) 
Xhatminus <- matrix(1,(m^2)*p,1) #Prior coefficient estimates. "hat" represents estimation.
Q <- diag((m^2)*p)*10^(-3)
R<-diag(2)
Pmin = diag((m^2)*p) 
Xhatstorage <- matrix(0,(m^2)*p,sL) #Matrix to store all X(t) coefficients
Yhat<-matrix(0,m,(sL-p)) #Matrix to store all VAR estimated data values
ch1<-1
ch2<-2

for(j in (p+1):(sL-p)){
  H <- diag(m) %x% t(vec(t(Y[c(ch1,ch2),(j-p):(j-1)])))
  K <- (Pmin %*% t(H)) %*% solve((H%*%Pmin%*%t(H) + R)) 
  Xhatplus <- F%*%(Xhatminus + K%*%(Y[c(ch1,ch2),j]-H%*%Xhatminus))
  E <- Y[c(ch1,ch2),j]-H%*%Xhatminus
  R <- 0.03*(E %*% t(E))
  P <- Pmin - K%*%H%*%Pmin
  Pplus <- (F%*% P %*% F) + Q
  Xhatminus <- Xhatplus
  Xhatstorage[,j]<- Xhatplus 
  Yhat[,j]<-H%*%Xhatplus
  Pmin <- Pplus
}

plot(1:(sL-p),Y[1,1:(sL-p)])  
lines(1:(sL-p),Yhat[1,1:(sL-p)],col="red") #See that VAR fits almost perfectly with data

plot(1:(sL-p),Y[2,1:(sL-p)])  
lines(1:(sL-p),Yhat[2,1:(sL-p)],col="red") #See that VAR fits almost perfectly with data

    plot(Xhatstorage[5,]) 
#Show coefficients from ch1->ch2 at t-2, visible influence at index 700-850, as expected from data simulation (from 700-800 values are t-2 from first time-series with added low variance noize)

    plot(Xhatstorage[6,]) 
#Show coefficients from ch1->ch2 at t-1, visible negative influence at index 300-700 as expected from data simulation (from 300-700 values are averages of first time-series multiplied by (-1))