diff --git a/tests/kalman/likelihood/compare_kalman_mex.m b/tests/kalman/likelihood/compare_kalman_mex.m
index 5655efad60e4f270dc1650d319539649841aebeb..1f9470ea0ea269e932bdc629badd56f6f82a0ced 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!');
-end
-
if measurement_error_flag==0
- HH = 0;
+ HH = [];
elseif measurement_error_flag==1
- HH = diag(H);
+ HH = diag(H);
elseif measurement_error_flag==2
- HH = H;
+ HH = H;
+end
+
+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;
-
-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)
+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;
+
+dprintf('Zflag = 0: Fortran Kalman filter is %5.2f times faster than its Matlab counterpart.', T_matlab_0/T_mex_0)
-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;
+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
+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;
-
-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)
+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 ((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;
+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_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;
+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
+
+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 92d3450e06da8c3f30fd49be25b6a833e06db666..8f0a2232709c96b9b14500f2d3c8d7c33e3cc00f 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 e039340963321995fde76737656186e652f3a19b..e1a1375ae06230fa8b3e71c1f11acc87c37359dd 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);