问题描述
我想使用MARSS（或R中的其他软件包）估算（MLE）该模型
x_t=x_{t-1}+w_t,with w_t ~ N(0,q)
y_t= d1_t + \alpha d2_t + \beta (d3_t -x_{t-1}) + v_t,with v_t ~ N(0,6*q)
其中第一行是转移方程，第二行是观察方程。 我设法以MARSS（R-package）接受的形式编写了它，如下所示：
[x1_t,x2_{t-1}]= [1,0;1,0][x1_{t-1},x2_{t-2}]+[w1_t,w2_t],with w1_t ~ N(0,q) and w2_t ~ N(0,0)
y_t= D d_t+Z x_t,6*q)
where
x_t=[x1_t,x2_{t-1}]
D=[1,\alpha,\beta]
Z=[0,\beta]
d_t=[d1_t,d2_t,d3_t]
问题是我无法使约束正常工作。当我运行此系统时，R会将Z矩阵中的\ beta与D矩阵中的\ beta分开考虑。我在互联网上看到的所有示例均显示仅使用Z矩阵（或仅使用D）的线性限制。我想做倍数的方差也会发生同样的情况。 有人可以帮我吗？
这是一个玩具数据：
B <- matrix(list(1,1,0),2,byrow=TRUE)
U <- matrix(0,1)
C <- matrix(0,1)
G <- matrix(list(1,byrow=TRUE)
Q <- matrix(list('d',byrow=TRUE)
Z <- matrix(list(0,'b'),2)
A <- matrix(0)
D <- matrix(list(1,'a',3)
H <- matrix(1)
R=matrix(list('6*d'))
dt<-matrix(rnorm(300),3,100)
y<-rnorm(100)
x0=matrix(list(0.094,0.094),1)
V0=matrix(list(0.001,0.001),2)
model.list = list(B=B,U=U,C=C,Q=Q,Z=Z,A=A,D=D,d=dt,H=H,R=R,x0=x0,V0=V0)
kemfit = MARSS(y,model=model.list,control=list(maxit=100,conv.test.slope.tol=0.1,abstol=0.1),method='kem')
解决方法
MARSS中的EM算法仅允许在相同矩阵内进行约束（例如设置值相等）。在A＆D或U＆C上设置约束很容易，但是在D＆Z或R＆Q上设置约束要求以一种奇怪的方式重写模型，其中您的协变量（dt）显示为虚拟状态（x）。所以你不想那样做。
您只需编写一个函数即可返回状态空间模型的负对数似然率，然后使用optim()将其最小化。我会使用SSCustom()函数对KFAS软件包进行此操作，因为那样会很快。但是，这里是如何使用MARSS软件包执行此操作的目的，只是向您展示了这一概念。作为MARSS的作者，我可以立即将其写下来，而使用KFAS软件包（我也使用过），我需要查找如何做协变量。
# Set up the parts that don't change
dt<-matrix(rnorm(300),3,100)
y<-rnorm(100)
x0=matrix(list(0.094,0.094),2,1)
V0=matrix(list(0.001,0.001),2)
B <- matrix(list(1,1,0),byrow=TRUE)
U <- A <- "zero"
# Put the parameters you will estimate into a vector
pars <- c(a=0.1624,b=-0.1,d=sqrt(0.2))
# Write a function to return the negative log-likelihood
negloglik <- function(pars){
Q <- matrix(list(pars["d"]^2,byrow=TRUE)
Z <- matrix(list(0,pars["b"]),2)
D <- matrix(list(1,pars["a"],3)
R <- matrix(6*pars["d"]^2)
model.list = list(B=B,U=U,Q=Q,Z=Z,A=A,D=D,d=dt,R=R,x0=x0,V0=V0)
-1*MARSS(y,model=model.list,control=list(maxit=100,conv.test.slope.tol=0.1,abstol=0.1),method='kem',silent=TRUE)$logLik
}
optim(pars,negloglik,method="BFGS")
在这里使用MARSS()函数获取logLik有点愚蠢，因为这是一个拟合函数，但是在固定了所有参数的情况下，它只会返回不拟合的logLik。
如果要查看KFAS模型的外观，可以执行以下操作：
kfas.model <- MARSSkfas(kemfit,return.kfas.model=TRUE,return.lag.one=FALSE)$kfas.model
然后
library(KFAS)
logLik(kfas.model)
将使您获得对数可能性。但是协变量如何进入KFAS模型有点不直观。它们在kfas.model$Z元素中显示为随时间变化的Z。我确信KFAS程序包具有一些辅助函数来构造带有协变量的模型。我总是从矩阵（没有辅助函数）构造KFAS模型，所以我对它们不熟悉，但是我知道它们存在。