Fminunc - complex standard errrors from taking inverse of Hessian
16 views (last 30 days)
Show older comments
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
0 Comments
Answers (2)
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.
Matt J
on 20 May 2019
You can also try this finite differencer to see if it gives you better Hessian calculation
https://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation
Using a finite difference approximation of the Hessian is not ideal, but this at least should be less crude than what BFGS computes.
0 Comments
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!