The example is based on an example in Shumway and Stoffer's: "Time Series Analysis and it's Applications with R Examples". In the original example phi, cq, and cr were scalar so the authors could use fdHess without any issues (see the hashed out version of the code).
para=list(phi, cq, cr)
Linn=function(para){# to evaluate likelhood at estimates
#kf=Kfilter0(num,y,1,mu0,Sigma0,para[1],para[2],para[3])
kf=Kfilter0(num,y,A=h,mu0,Sigma0,para[[1]],para[[2]],para[[3]])
return(kf$like)}
emhess=fdHess(para, function(para) Linn(para))
SE=sqrt(diag(solve(emhess$Hessian)))
I would like to generalize the code so that it can be applied to multivariate time series models. So in the code shown phi, cq, and cr are n*n arrays.
Is there a package that can calculate the Hessian for a scalar valued function with matrix arguments?
The closest match I can find is this (I also looked at nlme and numDeriv):
calculating the Gradient and the Hessian in R
In this case all the arguments are passed as a vector so the function being called has to be modified so that it can take the list of arguments and reconstruct the required matrices.
Is there a method that would allow me to calculate the Hessian for a scalar valued function with matrix arguments without changing the function being called? Seems this would be such a common problem that there would be an off the shelf answer but I haven't been able to find one.
Baz