replaced expressions X == 0 by more efficient and robust isequal(X,0) when X is a matrix
porting earlier changes from DsgeLikelihood.m to DsgeLikelihood_hh. cleaning codetime-shift
Michel Juillard
parent
f7477e23e5
commit
b247c329d8
|
|
@ -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...
|
||||
|
|
|
|||
|
|
@ -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...
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Reference in New Issue