diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m
index 6a4d2505c31d246541d27a4fe2bbdd3098410021..4468b411a61ccc4bef0a5d0393c2baf902c74b11 100644
--- a/matlab/DsgeLikelihood.m
+++ b/matlab/DsgeLikelihood.m
@@ -206,7 +206,7 @@ if (kalman_algo==1)% Multivariate Kalman Filter
end
if (kalman_algo==2)% Univariate Kalman Filter
no_correlation_flag = 1;
- if length(H)==1 & H == 0
+ if isequal(H,0)
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
@@ -237,7 +237,7 @@ if (kalman_algo==3)% Multivariate Diffuse Kalman Filter
end
if (kalman_algo==4)% Univariate Diffuse Kalman Filter
no_correlation_flag = 1;
- if length(H)==1 & H == 0
+ if isequal(H,0)
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
diff --git a/matlab/DsgeLikelihood_hh.m b/matlab/DsgeLikelihood_hh.m
index 8fc6629ea0da3b90722f61dcefe7aeb690bd53cc..e68370f114f90406cd64003bcc26cb7d084a050e 100644
--- a/matlab/DsgeLikelihood_hh.m
+++ b/matlab/DsgeLikelihood_hh.m
@@ -181,65 +181,7 @@ elseif options_.lik_init == 3 % Diffuse Kalman filter
if kalman_algo ~= 4
kalman_algo = 3;
end
- [QT,ST] = schur(T);
- e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
- [QT,ST] = ordschur(QT,ST,e1);
- k = find(abs(ordeig(ST)) > 2-options_.qz_criterium);
- nk = length(k);
- nk1 = nk+1;
- Pinf = zeros(np,np);
- Pinf(1:nk,1:nk) = eye(nk);
- Pstar = zeros(np,np);
- B = QT'*R*Q*R'*QT;
- for i=np:-1:nk+2
- if ST(i,i-1) == 0
- if i == np
- c = zeros(np-nk,1);
- else
- c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
- ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
- end
- q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
- Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
- Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
- else
- if i == np
- c = zeros(np-nk,1);
- c1 = zeros(np-nk,1);
- else
- c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
- ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
- ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
- c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
- ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
- ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
- end
- q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
- -ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
- z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
- Pstar(nk1:i,i) = z(1:(i-nk));
- Pstar(nk1:i,i-1) = z(i-nk+1:end);
- Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
- Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
- i = i - 1;
- end
- end
- if i == nk+2
- c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
- Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
- end
- Z = QT(mf,:);
- R1 = QT'*R;
- [QQ,RR,EE] = qr(Z*ST(:,1:nk),0);
- k = find(abs(diag([RR; zeros(nk-size(Z,1),size(RR,2))])) < 1e-8);
- if length(k) > 0
- k1 = EE(:,k);
- dd =ones(nk,1);
- dd(k1) = zeros(length(k1),1);
- Pinf(1:nk,1:nk) = diag(dd);
- end
-end
-if kalman_algo == 2
+ [Z,ST,R1,QT,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,options_.qz_criterium);
end
kalman_tol = options_.kalman_tol;
riccati_tol = options_.riccati_tol;
@@ -262,7 +204,7 @@ if (kalman_algo==1)% Multivariate Kalman Filter
end
if (kalman_algo==2)% Univariate Kalman Filter
no_correlation_flag = 1;
- if length(H)==1 & H == 0
+ if isequal(H,0)
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
@@ -293,7 +235,7 @@ if (kalman_algo==3)% Multivariate Diffuse Kalman Filter
end
if (kalman_algo==4)% Univariate Diffuse Kalman Filter
no_correlation_flag = 1;
- if length(H)==1 & H == 0
+ if isequal(H,0)
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 2c7a57a1f249bdc88a0f11f5d57a2aa54e2186ed..82fbf4db3b3090445d7c0056a3ba7a00193a6210 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -133,7 +133,7 @@ data1 = Y-trend;
% -----------------------------------------------------------------------------
% 4. Kalman smoother
% -----------------------------------------------------------------------------
-if size(H,1) == 1 && H == 0
+if isequal(H,0)
H = zeros(nobs,nobs);
end
@@ -156,7 +156,7 @@ if kalman_algo == 1 || kalman_algo == 3
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH1_Z(ST, ...
Z,R1,Q,H,Pinf,Pstar, ...
data1,nobs,np,smpl,kalman_tol,riccati_tol,data_index);
- if all(alphahat(:)==0)
+ if isequal(alphahat,0)
if kalman_algo == 1
kalman_algo = 2;
elseif kalman_algo == 3
@@ -214,147 +214,7 @@ if estim_params_.ncn && (kalman_algo == 2 || kalman_algo == 4)
if ~isempty(decomp)
decomp = decomp(:,k,:,:);
end
- if ~isempte(P)
+ if ~isempty(P)
P = P(k,k,:);
end
end
-
-% $$$ if any(any(H ~= 0)) % should be replaced by a flag
-% $$$ if kalman_algo == 1
-% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = ...
-% $$$ kalman_smoother(T,R,Q,H,Pstar,data1,start,mf,kalman_tol,riccati_tol);
-% $$$ if all(alphahat(:)==0)
-% $$$ kalman_algo = 2;
-% $$$ if ~estim_params_.ncn
-% $$$ [alphahat,epsilonhat,etahat,ahat,aK] = ...
-% $$$ DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
-% $$$ else
-% $$$ [alphahat,epsilonhat,etahat,ahat,aK] = ...
-% $$$ DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf,Pstar,Y,trend, ...
-% $$$ nobs,np,smpl,mf);
-% $$$ end
-% $$$ end
-% $$$ elseif kalman_algo == 2
-% $$$ if ~estim_params_.ncn
-% $$$ [alphahat,epsilonhat,etahat,ahat,aK] = ...
-% $$$ DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
-% $$$ else
-% $$$ [alphahat,epsilonhat,etahat,ahat,aK] = ...
-% $$$ DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf,Pstar,Y,trend, ...
-% $$$ nobs,np,smpl,mf);
-% $$$ end
-% $$$ elseif kalman_algo == 3 | kalman_algo == 4
-% $$$ data1 = Y - trend;
-% $$$ if kalman_algo == 3
-% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
-% $$$ DiffuseKalmanSmootherH1_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,nobs,np,smpl);
-% $$$ if all(alphahat(:)==0)
-% $$$ kalman_algo = 4;
-% $$$ if ~estim_params_.ncn
-% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
-% $$$ DiffuseKalmanSmootherH3_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,nobs,np,smpl);
-% $$$ else
-% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
-% $$$ DiffuseKalmanSmootherH3corr_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1, ...
-% $$$ nobs,np,smpl);
-% $$$ end
-% $$$ end
-% $$$ else
-% $$$ if ~estim_params_.ncn
-% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
-% $$$ DiffuseKalmanSmootherH3_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1, ...
-% $$$ nobs,np,smpl);
-% $$$ else
-% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,d,decomp] = ...
-% $$$ DiffuseKalmanSmootherH3corr_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1, ...
-% $$$ nobs,np,smpl);
-% $$$ end
-% $$$ end
-% $$$ alphahat = QT*alphahat;
-% $$$ ahat = QT*ahat;
-% $$$ nk = options_.nk;
-% $$$ for jnk=1:nk
-% $$$ aK(jnk,:,:) = QT*squeeze(aK(jnk,:,:));
-% $$$ for i=1:size(PK,4)
-% $$$ PK(jnk,:,:,i) = QT*squeeze(PK(jnk,:,:,i))*QT';
-% $$$ end
-% $$$ for i=1:size(decomp,4)
-% $$$ decomp(jnk,:,:,i) = QT*squeeze(decomp(jnk,:,:,i));
-% $$$ end
-% $$$ end
-% $$$ for i=1:size(P,4)
-% $$$ P(:,:,i) = QT*squeeze(P(:,:,i))*QT';
-% $$$ end
-% $$$ end
-% $$$ else
-% $$$ H = 0;
-% $$$ if kalman_algo == 1
-% $$$ if missing_value
-% $$$ [alphahat,etahat,ahat,aK] = missing_DiffuseKalmanSmoother1(T,R,Q, ...
-% $$$ Pinf,Pstar,Y,trend,nobs,np,smpl,mf,data_index);
-% $$$ else
-% $$$ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = ...
-% $$$ kalman_smoother(T,R,Q,H,Pstar,data1,start,mf,kalman_tol,riccati_tol);
-% $$$ end
-% $$$ if all(alphahat(:)==0)
-% $$$ kalman_algo = 2;
-% $$$ end
-% $$$ end
-% $$$ if kalman_algo == 2
-% $$$ if missing_value
-% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother3(T,R,Q, ...
-% $$$ Pinf,Pstar,Y,trend,nobs,np,smpl,mf,data_index);
-% $$$ else
-% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3(T,R,Q, ...
-% $$$ Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
-% $$$ end
-% $$$ end
-% $$$ if kalman_algo == 3
-% $$$ data1 = Y - trend;
-% $$$ if missing_value
-% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother1_Z(ST, ...
-% $$$ Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl,data_index);
-% $$$ else
-% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother1_Z(ST, ...
-% $$$ Z,R1,Q,Pinf,Pstar, ...
-% $$$ data1,nobs,np,smpl);
-% $$$ end
-% $$$ if all(alphahat(:)==0)
-% $$$ kalman_algo = 4;
-% $$$ end
-% $$$ end
-% $$$ if kalman_algo == 4
-% $$$ data1 = Y - trend;
-% $$$ if missing_value
-% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmoother3_Z(ST, ...
-% $$$ Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl,data_index);
-% $$$ else
-% $$$ [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ...
-% $$$ Z,R1,Q,Pinf,Pstar, ...
-% $$$ data1,nobs,np,smpl);
-% $$$ end
-% $$$ end
-% $$$ if kalman_algo == 3 | kalman_algo == 4
-% $$$ alphahat = QT*alphahat;
-% $$$ ahat = QT*ahat;
-% $$$ nk = options_.nk;
-% $$$ % $$$ if M_.exo_nbr<2 % Fix the crash of Dynare when the estimated model has only one structural shock (problem with
-% $$$ % $$$ % the squeeze function, that does not affect 2D arrays).
-% $$$ % $$$ size_decomp = 0;
-% $$$ % $$$ else
-% $$$ % $$$ size_decomp = size(decomp,4);
-% $$$ % $$$ end
-% $$$ for jnk=1:nk
-% $$$ aK(jnk,:,:) = QT*squeeze(aK(jnk,:,:));
-% $$$ for i=1:size(PK,4)
-% $$$ PK(jnk,:,:,i) = QT*dynare_squeeze(PK(jnk,:,:,i))*QT';
-% $$$ end
-% $$$ for i=1:size(decomp,4)
-% $$$ decomp(jnk,:,:,i) = QT*dynare_squeeze(decomp(jnk,:,:,i));
-% $$$ end
-% $$$ end
-% $$$ for i=1:size(P,4)
-% $$$ P(:,:,i) = QT*dynare_squeeze(P(:,:,i))*QT';
-% $$$ end
-% $$$ end
-% $$$ end
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index cc2e8a3b6d1eeb4fd9ec4bee9f8558a0e3e83bbb..9898651832cd855e9217c3038e10035c28c02c5a 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -220,7 +220,7 @@ else
k3 = (1:M_.endo_nbr)';
end
bayestopt_.smoother_var_list = union(k2,k3);
-[junk,bayestopt_.smoother_saved_var_list] = intersect(k3,bayestopt_.smoother_var_list);
+[junk,bayestopt_.smoother_saved_var_list] = intersect(k3,bayestopt_.smoother_var_list(:));
[junk,ic] = intersect(bayestopt_.smoother_var_list,nstatic+(1:npred)');
bayestopt_.smoother_restrict_columns = ic;
[junk,bayestopt_.smoother_mf] = ismember(var_obs_index, ...
@@ -306,7 +306,7 @@ options_ = set_default_option(options_,'mh_nblck',2);
options_ = set_default_option(options_,'nodiagnostic',0);
% load mode file is necessary
-if length(options_.mode_file) > 0 && ~options_.mh_posterior_mode_estimation
+if ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
load(options_.mode_file);
end
@@ -366,7 +366,7 @@ missing_value = ~(number_of_observations == gend*n_varobs);
initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations);
-if isnumeric(options_.mode_compute) && options_.mode_compute == 0 && length(options_.mode_file) == 0 && options_.mh_posterior_mode_estimation==0
+if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
if options_.smoother == 1
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,gend,data,data_index,missing_value);
oo_.Smoother.SteadyState = ys;
@@ -401,14 +401,14 @@ if isnumeric(options_.mode_compute) && options_.mode_compute == 0 && length(opti
return;
end
-if options_.mode_compute==6
+if isequal(options_.mode_compute,6)
% Erase previously computed optimal mh scale parameter.
delete([M_.fname '_optimal_mh_scale_parameter.mat'])
end
%% Estimation of the posterior mode or likelihood mode
-if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation
+if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
if ~options_.dsge_var
fh=str2func('DsgeLikelihood');
else
@@ -645,8 +645,7 @@ if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation
' option is unknown!'])
end
end
- % if options_.mode_compute ~= 5
- if options_.mode_compute ~= 6
+ if ~isequal(options_.mode_compute,6)
if options_.cova_compute == 1
if ~options_.dsge_var
hh = reshape(hessian('DsgeLikelihood',xparam1, ...
@@ -665,7 +664,6 @@ if any(options_.mode_compute ~= 0) && ~options_.mh_posterior_mode_estimation
else
save([M_.fname '_mode.mat'],'xparam1','parameter_names');
end
- % end
end
if options_.cova_compute == 0
diff --git a/matlab/set_prior.m b/matlab/set_prior.m
index 650ba9936b9262636db56f5efddbae2204dbaf80..b27af0248bb27004a3036d9ba4045c61af6ad16c 100644
--- a/matlab/set_prior.m
+++ b/matlab/set_prior.m
@@ -74,7 +74,7 @@ if nvx
bayestopt_.name = cellstr(M_.exo_names(estim_params_.var_exo(:,1),:));
end
if nvn
- if M_.H == 0
+ if isequal(M_.H,0)
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
end
@@ -115,7 +115,7 @@ if ncx
end
end
if ncn
- if M_.H == 0
+ if isequal(M_.H,0)
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
end