From b863c309bd5c961aa355fb4f92f39a8070d426ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Argos=29?= Date: Thu, 9 Nov 2023 17:06:26 +0100 Subject: [PATCH] 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). --- tests/kalman/likelihood/compare_kalman_mex.m | 134 ++++++++++-------- .../likelihood/simul_state_space_model.m | 17 +-- tests/kalman/likelihood/test_kalman_mex.m | 7 +- 3 files changed, 88 insertions(+), 70 deletions(-) diff --git a/tests/kalman/likelihood/compare_kalman_mex.m b/tests/kalman/likelihood/compare_kalman_mex.m index 5655efad6..1f9470ea0 100644 --- a/tests/kalman/likelihood/compare_kalman_mex.m +++ b/tests/kalman/likelihood/compare_kalman_mex.m @@ -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_matlab1e-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_matlab1e-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) diff --git a/tests/kalman/likelihood/simul_state_space_model.m b/tests/kalman/likelihood/simul_state_space_model.m index 92d3450e0..8f0a22327 100644 --- a/tests/kalman/likelihood/simul_state_space_model.m +++ b/tests/kalman/likelihood/simul_state_space_model.m @@ -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 \ No newline at end of file + end diff --git a/tests/kalman/likelihood/test_kalman_mex.m b/tests/kalman/likelihood/test_kalman_mex.m index e03934096..e1a1375ae 100644 --- a/tests/kalman/likelihood/test_kalman_mex.m +++ b/tests/kalman/likelihood/test_kalman_mex.m @@ -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);