Improve tests and timings for the Kalman mex.
- Ensure that we always use the same DGP (e.g. same transition matrix). - Call the mex more than once with different samples. - Ensure that the size of the state space model is the same in tests 1, 2 and 3. - Fix the seed (same samples across 1, 2 and 3 except for the additive noise in 2 and 3 on the observed variables).kalman_mex
parent
9c61422a75
commit
b863c309bd
GPG Key ID:
A6D44CB9C64CE77B
|
|
@ -1,31 +1,39 @@
|
|||
function [flag] = compare_kalman_mex(experience)
|
||||
|
||||
pp = experience.Number0fObservedVariables;
|
||||
mm = experience.SizeOfTheStateVector;
|
||||
rr = experience.NumberOfStructuralShocks;
|
||||
|
||||
rng(experience.Seed);
|
||||
|
||||
N = experience.NumberOfSimulations;
|
||||
|
||||
pp = experience.Number0fObservedVariables*experience.Scale;
|
||||
mm = experience.SizeOfTheStateVector*experience.Scale;
|
||||
rr = experience.NumberOfStructuralShocks*experience.Scale;
|
||||
measurement_error_flag = experience.MeasurementErrors;
|
||||
gend = experience.NumberOfPeriods;
|
||||
|
||||
%% SET VARIOUS PARAMETERS
|
||||
%% SET VARIOUS PARAMETERS
|
||||
kalman_tol = 1e-12;
|
||||
riccati_tol =1e-9;
|
||||
|
||||
%% SET THE STATE SPACE MODEL:
|
||||
|
||||
% I randomly choose the mm eigenvalues of the transition matrix.
|
||||
TransitionEigenvalues = rand(mm,1)*2-1;
|
||||
TransitionEigenvalues = rand(mm, 1)*2-1;
|
||||
|
||||
% I randomly choose the mm eigenvectors of the transition matrix
|
||||
tmp = rand(mm,mm*100);
|
||||
TransitionEigenvectors = tmp*tmp'/(mm*100);
|
||||
TransitionEigenvectors = rand(mm,mm);
|
||||
TransitionEigenvectors = rand(mm, mm);
|
||||
|
||||
% I build the transition matrix
|
||||
T = TransitionEigenvectors*diag(TransitionEigenvalues)/TransitionEigenvectors;
|
||||
|
||||
% I randomly choose matrix R
|
||||
R = randn(mm,rr);
|
||||
R = randn(mm, rr);
|
||||
|
||||
% I randomly choose the covariance matrix of the structural innovations
|
||||
E = randn(rr,20*rr);
|
||||
E = randn(rr, 20*rr);
|
||||
Q = E*transpose(E)/(20*rr);
|
||||
% If needed I randomly choose the covariance matrix of the measurement errors
|
||||
|
||||
% If needed I randomly choose the covariance matrix of the measurement errors
|
||||
if measurement_error_flag == 0
|
||||
H = zeros(pp,1);
|
||||
elseif measurement_error_flag == 1
|
||||
|
|
@ -33,82 +41,92 @@ elseif measurement_error_flag == 1
|
|||
elseif measurement_error_flag == 2
|
||||
E = randn(pp,20*pp);
|
||||
H = E*transpose(E)/(20*pp);
|
||||
H = (.1*eye(pp))*H*(.1*eye(pp));
|
||||
else
|
||||
disp('compare_kalman_mex: unknown option!')
|
||||
end
|
||||
% Set the selection vector (mf)
|
||||
|
||||
% Set the selection vector (mf)
|
||||
MF = transpose(randperm(mm));
|
||||
mf = MF(1:pp);
|
||||
|
||||
% Compute ergodic variance of the state equation
|
||||
P = lyapunov_symm(T,R*Q*R',riccati_tol,1.000001, riccati_tol);
|
||||
|
||||
%% BUILD DATA SET (zero mean):
|
||||
a = zeros(mm,1);
|
||||
if measurement_error_flag == 0
|
||||
Y = simul_state_space_model(T,R,Q,mf,gend);
|
||||
elseif measurement_error_flag == 1
|
||||
H = rand(pp,1);
|
||||
Y = simul_state_space_model(T,R,Q,mf,gend,diag(H));
|
||||
elseif measurement_error_flag == 2
|
||||
E = randn(pp,20*pp);
|
||||
H = E*transpose(E)/(20*pp);
|
||||
Y = simul_state_space_model(T,R,Q,mf,gend,H);
|
||||
else
|
||||
disp('compare_kalman_mex: unknown option!');
|
||||
if measurement_error_flag==0
|
||||
HH = [];
|
||||
elseif measurement_error_flag==1
|
||||
HH = diag(H);
|
||||
elseif measurement_error_flag==2
|
||||
HH = H;
|
||||
end
|
||||
|
||||
if measurement_error_flag==0
|
||||
HH = 0;
|
||||
elseif measurement_error_flag==1
|
||||
HH = diag(H);
|
||||
elseif measurement_error_flag==2
|
||||
HH = H;
|
||||
rng(experience.Seed*1938);
|
||||
|
||||
% Build datasets
|
||||
Y = simul_state_space_model(T,R,Q,mf,gend*N, HH);
|
||||
|
||||
if isempty(HH)
|
||||
HH = 0;
|
||||
end
|
||||
|
||||
%
|
||||
% Evaluate likelihoods
|
||||
%
|
||||
|
||||
LIK_matlab_0 = NaN(N,1);
|
||||
LIK_mex_0 = NaN(N,1);
|
||||
|
||||
flag = 0;
|
||||
Zflag = 0;
|
||||
tic;
|
||||
[LIK_matlab,lik_matlab] = kalman_filter(Y,1,gend,a,P,kalman_tol,riccati_tol,0,0,T,Q,R,HH,mf,mm,pp,rr,Zflag,0,0);
|
||||
T_matlab = toc;
|
||||
for i=1:N
|
||||
[LIK_matlab_0(i), ~] = kalman_filter(Y(:,(i-1)*gend+1:i*gend),1,gend,zeros(mm,1),P,kalman_tol,riccati_tol,0,0,T,Q,R,HH,mf,mm,pp,rr,Zflag,0,0);
|
||||
end
|
||||
T_matlab_0 = toc;
|
||||
|
||||
tic;
|
||||
[LIK_mex,lik_mex] = kalman_filter_mex(Y,a,P,kalman_tol,riccati_tol,T,Q,R,mf,Zflag,HH);
|
||||
T_mex = toc;
|
||||
for i=1:N
|
||||
[LIK_mex_0(i), ~] = kalman_filter_mex(Y(:,(i-1)*gend+1:i*gend),zeros(mm,1),P,kalman_tol,riccati_tol,T,Q,R,mf,Zflag,HH);
|
||||
end
|
||||
T_mex_0 = toc;
|
||||
|
||||
if T_matlab<T_mex
|
||||
dprintf('Zflag = 0: Matlab Kalman filter is %5.2f times faster than its Fortran counterpart.', T_mex/T_matlab)
|
||||
else
|
||||
dprintf('Zflag = 0: Fortran Kalman filter is %5.2f times faster than its Matlab counterpart.', T_matlab/T_mex)
|
||||
dprintf('Zflag = 0: Fortran Kalman filter is %5.2f times faster than its Matlab counterpart.', T_matlab_0/T_mex_0)
|
||||
|
||||
if max(abs((LIK_mex_0-LIK_matlab_0)./LIK_matlab_0))>1e-6
|
||||
dprintf("Zflag = 0: discrepancy between Matlab and Fortran Kalman filter results!")
|
||||
flag = 1;
|
||||
end
|
||||
|
||||
if ((abs((LIK_matlab - LIK_mex)/LIK_matlab) > 1e-6) || (max(abs((lik_matlab-lik_mex)./lik_matlab)) > 1e-6))
|
||||
dprintf("Zflag = 0: discrepancy between Matlab and Fortran Kalman filter results!")
|
||||
flag = 1;
|
||||
end
|
||||
LIK_matlab_1 = NaN(N,1);
|
||||
LIK_mex_1 = NaN(N,1);
|
||||
|
||||
Zflag = 1;
|
||||
Z = eye(mm);
|
||||
Z = Z(mf,:);
|
||||
tic;
|
||||
[LIK_matlab_z,lik_matlab_z] = kalman_filter(Y,1,gend,a,P,kalman_tol,riccati_tol,0,0,T,Q,R,HH,Z,mm,pp,rr,Zflag,0,0);
|
||||
T_matlab = toc;
|
||||
for i=1:N
|
||||
[LIK_matlab_1(i), ~] = kalman_filter(Y(:,(i-1)*gend+1:i*gend),1,gend,zeros(mm,1),P,kalman_tol,riccati_tol,0,0,T,Q,R,HH,Z,mm,pp,rr,Zflag,0,0);
|
||||
end
|
||||
T_matlab_1 = toc;
|
||||
|
||||
tic;
|
||||
[LIK_mex_z,lik_mex_z] = kalman_filter_mex(Y,a,P,kalman_tol,riccati_tol,T,Q,R,Z,Zflag,HH);
|
||||
T_mex = toc;
|
||||
for i=1:N
|
||||
[LIK_mex_1(i), ~] = kalman_filter_mex(Y(:,(i-1)*gend+1:i*gend),zeros(mm,1),P,kalman_tol,riccati_tol,T,Q,R,Z,Zflag,HH);
|
||||
end
|
||||
T_mex_1 = toc;
|
||||
|
||||
if T_matlab<T_mex
|
||||
dprintf('Zflag = 1: Matlab Kalman filter is %5.2f times faster than its Fortran counterpart.', T_mex/T_matlab)
|
||||
else
|
||||
dprintf('Zflag = 1: Fortran Kalman filter is %5.2f times faster than its Matlab counterpart.', T_matlab/T_mex)
|
||||
dprintf('Zflag = 1: Fortran Kalman filter is %5.2f times faster than its Matlab counterpart.', T_matlab_1/T_mex_1)
|
||||
|
||||
if max(abs((LIK_mex_1-LIK_matlab_1)./LIK_matlab_1))>1e-6
|
||||
dprintf("Zflag = 1: discrepancy between Matlab and Fortran Kalman filter results!")
|
||||
flag = 1;
|
||||
end
|
||||
|
||||
if ((abs((LIK_matlab_z - LIK_mex_z)/LIK_matlab_z) > 1e-6) || (max(abs((lik_matlab_z-lik_mex_z)./lik_matlab_z)) > 1e-6))
|
||||
dprintf("Zflag = 1: discrepancy between Matlab and Fortran Kalman filter results!")
|
||||
flag = 1;
|
||||
if max(abs((LIK_mex_0 - LIK_mex_1)./LIK_mex_1))>1e-6
|
||||
dprintf("Zflag = 1: discrepancy between results with and without the Zflag!")
|
||||
flag = 1;
|
||||
end
|
||||
|
||||
if ((abs((LIK_mex - LIK_mex_z)/LIK_mex) > 1e-6) || (max(abs((lik_mex-lik_mex_z)./lik_mex)) > 1e-6))
|
||||
dprintf("Zflag = 1: discrepancy between results with and without the Zflag!")
|
||||
flag = 1;
|
||||
end
|
||||
dprintf('Cost of Zflag==1 is %5.2f%% (matlab)', 100*(T_matlab_1-T_matlab_0)/T_matlab_0)
|
||||
dprintf('Cost of Zflag==1 is %5.2f%% (mex)', 100*(T_mex_1-T_mex_0)/T_mex_0)
|
||||
|
|
|
|||
|
|
@ -2,24 +2,25 @@ function observed_data = simul_state_space_model(T,R,Q,mf,nobs,H)
|
|||
pp = length(mf);
|
||||
mm = length(T);
|
||||
rr = length(Q);
|
||||
|
||||
|
||||
upper_cholesky_Q = chol(Q);
|
||||
if nargin>5
|
||||
if nargin>5 && ~isempty(H)
|
||||
upper_cholesky_H = chol(H);
|
||||
end
|
||||
|
||||
|
||||
state_data = zeros(mm,1);
|
||||
|
||||
if (nargin==5)
|
||||
|
||||
if (nargin==5 || isempty(H))
|
||||
for t = 1:nobs
|
||||
state_data = T*state_data + R* upper_cholesky_Q * randn(rr,1);
|
||||
observed_data(:,t) = state_data(mf);
|
||||
end
|
||||
elseif (nargin==6)
|
||||
elseif (nargin==6 && ~isempty(H))
|
||||
for t = 1:nobs
|
||||
state_data = T*state_data + R* upper_cholesky_Q * randn(rr,1);
|
||||
observed_data(:,t) = state_data(mf) + upper_cholesky_H * randn(pp,1);
|
||||
observed_data(:,t) = state_data(mf);
|
||||
end
|
||||
observed_data = observed_data + upper_cholesky_H * randn(pp,nobs);
|
||||
else
|
||||
error('simul_state_space_model:: I don''t understand what you want!!!')
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -23,6 +23,9 @@ Experience.SizeOfTheStateVector = 100;
|
|||
Experience.NumberOfStructuralShocks = 12;
|
||||
Experience.MeasurementErrors = 0;
|
||||
Experience.NumberOfPeriods = 300;
|
||||
Experience.Seed = 1;
|
||||
Experience.NumberOfSimulations = 100;
|
||||
Experience.Scale = 1;
|
||||
|
||||
try
|
||||
flag = compare_kalman_mex(Experience);
|
||||
|
|
@ -58,11 +61,7 @@ catch
|
|||
end
|
||||
|
||||
dprintf('Test 3: measurement error with general variance-covariance matrix')
|
||||
Experience.Number0fObservedVariables = 50;
|
||||
Experience.SizeOfTheStateVector = 70;
|
||||
Experience.NumberOfStructuralShocks = 50;
|
||||
Experience.MeasurementErrors = 2;
|
||||
Experience.NumberOfPeriods = 300;
|
||||
|
||||
try
|
||||
flag = compare_kalman_mex(Experience);
|
||||
|
|
|
|||
Loading…
Reference in New Issue