Fminunc - complex standard errrors from taking inverse of Hessian

16 views (last 30 days)
Does anyone know how to avoid having complex standard errors when using fminunc? The code convergences but at the minimum the standard errors are imaginary... I replicate the paper of De Jong (2000) who applies a kalman filter to esimate affine term structure models, this is the one-factor vasicek case. I hope someone knows how to solve this, I have provided my main script and the function used for optimization. Thank you!
Main code:
% ---- clear and load yields--------------------------------------------
clear
Y_load = xlsread('de_jong_1monthto10yeardata');
dt =1/12;
% Four maturity data set settings:
Y = [Y_load(:,5),Y_load(:,14),Y_load(:,26),Y_load(:,31)]'/100;
tau = [1/4 1 5 10];
tau = tau';
[nrow, ncol] = size(Y);
A_0 = 0.11;
kappa1 = 0.023; kappa2;
alpha1 = 0.00015; alpha2 = 0.00065;
psi1 = -2.75; psi2 = -16;
Z = [0.0006*ones(1,4),0.000001*ones(1,6)];
% Initialize parameters: take log to ensure positive variables
para0 = [A_0, kappa1,kappa2, alpha1, alpha2, psi1,psi2, Z];
para0(1:5) = log(para0(1:5));
para0(8:11) = log(para0(8:11));
% Optimize: [0.0604 0.0224 0.0002 -10.8003](complex hessian) 2lnL= -1.0035e+04
clearvars options
options = optimoptions('fminunc','Algorithm','quasi-newton','HessUpdate','bfgs','MaxFunctionEvaluations',1e+6,'MaxIterations',1e+6,'Display','iter');
[x,fval,exitflag,~,~,hessian] = fminunc(@(x) VAS_TWO_loglik_FUNC_NEW(x,Y,tau,ncol,nrow,dt),para0,options);
hes = hessian;
standarddeviatie(j,:) = sqrt(diag((1/ncol)*inv((hes))))';
Function:
function [sum_min_2_ln_L,F,P,predictedF,predictedP] = VAS_TWO_loglik_FUNC_NEW(x,Y,tau,ncol,nrow,dt) %calculate log likelihood
%% initialize the parameter for VAS model: De Jong check
% --- parameters -------------------------------------------------------
A_0 = exp(x(1));
kappa1 = exp(x(2)); kappa2 = exp(x(3));
alpha1 = exp(x(4)); alpha2 = exp(x(5));
psi1 = x(6); psi2 = x(7);
Z = [exp(x(8:11)),x(12:17)];
%% Create A and B: [KLOPT]
A = zeros(nrow,1);
B = zeros(nrow,2);
% one factor
AffineBeta1 = (ones(4,1) - exp(-kappa1*tau))/kappa1;
R_inf_1 = A_0 - (psi1*alpha1)/kappa1 - (alpha1/(2*kappa1^2));
Affinealpha1 = R_inf_1*(tau-AffineBeta1)+(alpha1/(4*kappa1))*AffineBeta1.^2;
A1 = Affinealpha1./tau;
B1 = AffineBeta1./tau;
% two factor
AffineBeta2 = (ones(4,1) - exp(-kappa2*tau))/kappa2;
R_inf_2 = A_0 - (psi2*alpha2)/kappa2 - (alpha2/(2*kappa2^2));
%R_inf = A_0 - (psi*sqrt(alpha_tilde))/kappa - (alpha_tilde/(2*kappa^2)); % dit is met sqrt want alpha_tilde = sigma^2 toch?
Affinealpha2 = R_inf_2*(tau-AffineBeta2)+(alpha2/(4*kappa2))*AffineBeta2.^2;
A2 = Affinealpha2./tau;
B2 = AffineBeta2./tau;
% Create A and B
A = A1+A2;
B(:,1) = B1;
B(:,2) = B2;
% --- Determine Phi -----------------------------------------------------
Phi = [exp(-kappa1*dt),0;0,exp(-kappa2*dt)];
% --- Determine Q -----------------------------------------------------
Q = [(alpha1/(2*kappa1))*(1-exp(-2*kappa1*dt)),0;0,(alpha2/(2*kappa2))*(1-exp(-2*kappa2*dt))];
% --- Determine H ------------------------------------------------------
L = [1,0,0,0;
Z(5), 1, 0, 0;
Z(6), Z(8), 1,0;
Z(7), Z(9),Z(10),1];
D = diag(Z(1:4));
H = L*D*L';
n=1;% number of factors
%% Kalman filter
% preallocate
min_2_ln_L = zeros(ncol,1); %[Check]
F = zeros(2,ncol);
P = zeros(2,n*ncol);
predictedF=zeros(2,ncol);
predictedP=zeros(2,n*ncol);
n=2;
for i = 1:ncol
% prediction:
if i==1
F(:,i) = [A_0;A_0];% F0 = E(Ft)
P(:,(n*i-1):(n*i)) = [alpha1/(2*kappa1),0;0,alpha2/(2*kappa2)];% P0 = var(Ft)
predictedF(:,i) = [A_0;A_0];
predictedP(:,(n*i-1):(n*i)) = [alpha1/(2*kappa1),0;0,alpha2/(2*kappa2)];
else
predictedF(:,i) = Phi * F(:,i-1);
predictedP(:,(n*i-1):(n*i)) = Phi * P(:,(n*i-3):(n*i-2)) * Phi' + Q;
end
% likelihood contributions [KLOPT]
u = Y(:,i) - A - B * predictedF(:,i); % error [4x1]
V = B*predictedP(:,(n*i-1):(n*i))*B'+ H;
if V<1e300; invv = pinv(V);
else invv=0;
end %in case kalman tries undefined values of pinv(V), substitute pinv(V) by 'invv' and uncomment
min_2_ln_L(i) = log(det(V)) + (transpose(u)*invv* u);
% updating [KLOPT]
K= predictedP(:,(n*i-1):(n*i))*B'*invv;
L = eye(n) - K*B;
F(:,i) = predictedF(:,i) + K*u;
P(:,(n*i-1):(n*i)) = L * predictedP(:,(n*i-1):(n*i));
P(:,(n*i-1):(n*i)) = (P(:,(n*i-1):(n*i))+ P(:,(n*i-1):(n*i))')/2;
end
sum_min_2_ln_L = sum(min_2_ln_L);
end

Answers (2)

Matt J
Matt J on 20 May 2019
Edited: Matt J on 20 May 2019
For the purpose of standard error computation, you should do a separate analytical computation of the Hessian at the final point. The Hessian returned by fminunc has lots of approximations in it, both due to the finite difference derivative calculations and because the BFGS algorithm may converge only semi-accurately to the actual Hessian.
You should also verify that the Hessian at the final point is positive definite with a good condition number cond(Hessian). That is the only situation when standard error approximation using Hessians is valid. Also, if the Hessian is not positive definite, it is possible that fminunc has converged to a saddle point, rather than to a local minimum.
  8 Comments
Tim Everaert
Tim Everaert on 20 May 2019
Matt! Thank you for your help!
I understand what you mean, so I manage to obtain the first function which just evaluates the LL function at the MLE estimates. I am not sure what you mean by the second? To be more clear, I have no clue how to write a function (Hessian(x)) in your example with regard to my function "VAS_TWO_loglik_FUNC_NEW"...
I am now also working on the finite differencer package you stated below, thanks for your help!
I would prefer the best option though, but still not see what you mean unfortunately...
Matt J
Matt J on 20 May 2019
Edited: Matt J on 20 May 2019
To be more clear, I have no clue how to write a function (Hessian(x)) in your example with regard to my function
Here's the link again that explains about Hessians

Sign in to comment.


Matt J
Matt J on 20 May 2019
You can also try this finite differencer to see if it gives you better Hessian calculation
Using a finite difference approximation of the Hessian is not ideal, but this at least should be less crude than what BFGS computes.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!