我正在使用卡尔曼滤波器来估计收益率曲线的各种动态和无套利纳尔逊-西格尔模型。我给出了一些起始值进行优化，算法收敛得很好。但是，当我想使用最优算法提供的 Hessian 计算标准误差时，由于方差协方差矩阵对角线上的非正值，我得到了 NaN。我认为这是因为我有一个具有许多局部最优值的高度非线性函数，但是对于我尝试的所有起始值，它都会不断发生。
我使用的函数与默认的 Nelder-Mead 算法一起optim。 我使用的命令是opt_para<-optim(par=par0, fn=Kalman_filter, y=y,
maturities=maturities,control=list(maxit=20000),hessian=TRUE)起始值以par0为单位给出，即
> par0
[1] 9.736930e-01 1.046646e+00 5.936238e-01 4.444669e-02 2.889251e-07 6.646960e+00 7.715964e-01 9.945551e-01 9.663361e-01
[10] 6.000000e-01 6.000000e-01 6.000000e-01 6.000000e-02 5.000000e-01 5.000000e-01 5.000000e-01 5.000000e-01
我得到的optim输出是
$
[1] 0.833208307 1.373442068 0.749313983 0.646577154 0.237102069 6.882644818 0.788775982 0.918378263 0.991982038
[10] 0.748509055 0.005115171 0.392213941 0.717186499 0.121525623 0.386227284
0.001970431 0.845279611$value
[1] 575.7886
$counts
function gradient
5225 NA
$convergence
[1] 0
$message
NULL
然后，我使用以下命令生成估计值的标准误差。
hessian<-opt_para$hessian
fish_info<-solve(hessian,tol=1e-100)
st_errors<- diag(sqrt(fish_info))
st_errors
我得到以下输出st_errors
[1] NaN NaN 2.9170315888 NaN NaN NaN 0.0294300357 0.0373614751 NaN
[10] 0.0785349634 0.0005656580 NaN 0.0470600219 0.0053255251 0.0408666177 0.0001561243 0.4540428740
NaN 在对角线上产生为负值，这在方差-协方差矩阵中是不可能的。但是，我怀疑这是由于优化过程不正确。
需要明确的是，我还包含要优化的功能。它是一个卡尔曼滤波器，内置了更新方程和一些限制。
Kalman_filter<-function(par, y, maturities){
b0<-c(par[1],par[2],par[3])
P0<-diag(c(par[4],par[5],par[6]))
Phi<-diag(c(par[7],par[8],par[9]))
mu<-c(par[10],par[11],par[12])
lambda<-par[13]
sigma11<-par[14]
sigma21<-par[15]
sigma22<-par[16]
sigma33<-par[17]
m=length(b0)
n=length(y[,1])
d<-length(y[1,])
sigma_eps<-sigma11*diag(d)
sigma_nu<-diag(c(sigma21^2,sigma22^2,sigma33^2))*(1/12)
colnames(sigma_nu)<-c("level","slope","Curvat")
X<-matrix(cbind(rep(1,length(maturities)), slope_factor(lambda,maturities), curv_factor(lambda,maturities)),ncol=3)
colnames(X)<-c("level","slope","Curvature")
bt<-matrix(NA, nrow=m, ncol=n+1)
Pt<-array(NA, dim=c(m,m,n+1))
btt<-matrix(NA, nrow=m,ncol=n+1)
Ptt<-array(NA, dim=c(m,m,n+1))
vt<-matrix(NA, nrow=d, ncol=n)
eigen_values<-eigen(Phi,only.values=TRUE)$valuesif(eigen_values[1]>=1||eigen_values[2]>=1||eigen_values[3]>=1){日志样=-70000000 }else{
c<- (diag(3) - Phi)%*% mu
loglike<-0i<-1btt[,1]<-b0Ptt[,,1]<-P0while(i< n+1){
bt[,i]<- c+ Phi%*% btt[,i]Pt[,,i] <- Phi%*% tcrossprod(Ptt[,,i],Phi) + sigma_nu
vt[,i]<- y[i,] - X%*% bt[,i]
ft<-X%*% tcrossprod(Pt[,,i], X) + sigma_eps
det_f<-det(ft)
if( is.nan(det_f) || is.na(det_f)|| is.infinite(det_f)){
loglike<- - 700000000
} else
{
if(det_f<0){
loglike <- - 700000000
} else
{
if (abs(det_f>1e-20)){
logdet_f<- log(det_f)
f_inv<- solve(ft, tol=1e-200)
Kt<- tcrossprod(Pt[,,i],X)%*% f_inv
btt[,i+1] <- bt[,i] + Kt%*% vt[,i]
Ptt[,,i+1] <- (diag(3) - Kt%*% X)%*% Pt[,,i]
loglike_contr<- -0.5*d*log(2*pi) - 0.5 * logdet_f - 0.5*
crossprod(vt[,i],f_inv)%*% vt[,i]
loglike<-loglike+loglike_contr
} else
{ loglike<- -700000}
}
}
i<-i+1
}
}
return(-loglike)
}
任何帮助将不胜感激。
我刚刚解决了这个问题，我再次使用唯一的输入参数(来自 optim 的似然估计(对似然函数进行了编程。在此之后，我使用了numDeriv包中的hessian函数。这为标准误差生成了可行的估计值。