#r #time-series #kalman-filter #state-space
#r #временные ряды #фильтр Калмана #пространство состояний
Вопрос:
В настоящее время я пытаюсь получить те же результаты оценки, что и в статье «Сезонные изменения температуры в Центральной Англии — Proietti amp; Hillebrand (2015)«
В этой модели они определяют следующую модель пространства состояний, которую я переписал, используя «стандартные» уравнения пространства состояний из Koopman amp; Durbin (2012). Синие переменные — это те, которые необходимо оценить, ожидаемые результаты приведены в первой статье (стр. 14).
Очень небольшое резюме для понимания модели: y_t — вектор размером 1xT с данными за месяц, X_t — вектор «выбора» размером 12×1, такой, что y_t допускает разные значения параметров за месяц, в котором он находится.
Я попытался сначала оценить модель, игнорируя матрицу X_t, которая изначально выдавала несколько ошибок… Например, первому элементу в Q не разрешалось быть равным нулю, что, на мой взгляд, странно.. у mu нет инноваций и, следовательно, нет различий для вектора инноваций, верно? Кроме того, даже при 1000 итерациях сходимость все еще не достигнута. Я думаю, что завтра утром я попытаюсь выполнить 10000 итераций, возможно, это поможет.
Однако главное, что меня ставит в тупик, это то, как разрешить X_t в этой модели пространства состояний.. это, в свою очередь, изменит модель на 60 параметров (поскольку каждый параметр представляет собой вектор размером 12×1). Любая помощь была бы прекрасной
library(KFAS)
yTOT = yTOT[yTOT$YR > 1771,]
yTOT = yTOT[yTOT$YR < 2014,]
y = t(as.matrix(yTOT$VAL))
B = matrix(list(0,"b1",0,0,1,0,0,0,"b2"),3,3, byrow=TRUE)
U = matrix(list("u1",0,0),3,1)
C = 0
c = 0
G = matrix(list(1,0,0, 0,1,0, 0,0,1),3,3, byrow=TRUE)
w = matrix(list(0,"w1","w2"),3,1)
Z = matrix(c(1,0,1),1,3)
A = matrix(0, 1,1)
D = 0
d = 0
H = 0
v = 0
R = matrix(0, 1,1)
Q = matrix(list(1,0,0, 0,1,0, 0,0,"q1"),3,3, byrow=TRUE)
x0 = matrix(0, 3,1)
V0 = matrix(0, 3,3)
model.gen=list(Z=Z,A=A,R=R,B=B,U=U,Q=Q,x0=x0,V0=V0,tinitx=0)
kemfit = MARSS(y, model=model.gen, control=list(maxit=1000))