diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index 729f5e43ca634c6e3024efc8fc9c6952bf072e57..ab678c33b889eacbac64cf566febe72edfda67d5 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -104,10 +104,6 @@ end
if options_.order>2 && options_.particle.pruning
error('Higher order nonlinear filters are not compatible with pruning option.')
end
-% Check the perturbation order (nonlinear filters with third order perturbation, or higher order, are not yet implemented).
-if options_.order>2 && ~isfield(options_,'options_ident') %For identification at order=3 we skip this check.
- error(['I cannot estimate a model with a ' int2str(options_.order) ' order approximation of the model!'])
-end
% analytical derivation is not yet available for kalman_filter_fast
if options_.analytic_derivation && options_.fast_kalman_filter
diff --git a/matlab/get_identification_jacobians.m b/matlab/get_identification_jacobians.m
index 8d66d438f9fbd57795b35dcedc864659c55e9df0..6481577654c401ca1c8bd9a1e88f52f3933544b7 100644
--- a/matlab/get_identification_jacobians.m
+++ b/matlab/get_identification_jacobians.m
@@ -220,13 +220,13 @@ options.options_ident.indpmodel = indpmodel;
options.options_ident.indpstderr = indpstderr;
options.options_ident.indpcorr = indpcorr;
oo.dr = pruned_state_space_system(M, options, oo.dr);
-E_y = oo.dr.pruned.E_y; dE_y = oo.dr.pruned.dE_y;
-A = oo.dr.pruned.A; dA = oo.dr.pruned.dA;
-B = oo.dr.pruned.B; dB = oo.dr.pruned.dB;
-C = oo.dr.pruned.C; dC = oo.dr.pruned.dC;
-D = oo.dr.pruned.D; dD = oo.dr.pruned.dD;
-c = oo.dr.pruned.c; dc = oo.dr.pruned.dc;
-d = oo.dr.pruned.d; dd = oo.dr.pruned.dd;
+E_y = oo.dr.pruned.E_y; dE_y = oo.dr.pruned.dE_y;
+A = oo.dr.pruned.A; dA = oo.dr.pruned.dA;
+B = oo.dr.pruned.B; dB = oo.dr.pruned.dB;
+C = oo.dr.pruned.C; dC = oo.dr.pruned.dC;
+D = oo.dr.pruned.D; dD = oo.dr.pruned.dD;
+c = oo.dr.pruned.c; dc = oo.dr.pruned.dc;
+d = oo.dr.pruned.d; dd = oo.dr.pruned.dd;
Varinov = oo.dr.pruned.Varinov; dVarinov = oo.dr.pruned.dVarinov;
Om_z = oo.dr.pruned.Om_z; dOm_z = oo.dr.pruned.dOm_z;
Om_y = oo.dr.pruned.Om_y; dOm_y = oo.dr.pruned.dOm_y;
@@ -259,13 +259,20 @@ if ~no_identification_moments
dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr);
dMOMENTS(1:obs_nbr,:) = dE_y; %add Jacobian of first moments of VAROBS variables
% Denote Ezz0 = E[zhat*zhat'], then the following Lyapunov equation defines the autocovariagram: Ezz0 -A*Ezz0*A' = B*Sig_xi*B' = Om_z
- Ezz0 = lyapunov_symm(A, B*Varinov*B', options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug);
- Eyy0 = C*Ezz0*C' + D*Varinov*D';
+ [Ezz0,u] = lyapunov_symm(A, Om_z, options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug);
+ stationary_vars = (1:length(indvobs))';
+ if ~isempty(u)
+ x = abs(C*u);
+ stationary_vars = find(all(x < options.Schur_vec_tol,2));
+ end
+ Eyy0 = NaN*ones(obs_nbr,obs_nbr);
+ Eyy0(stationary_vars,stationary_vars) = C(stationary_vars,:)*Ezz0*C(stationary_vars,:)' + Om_y(stationary_vars,stationary_vars);
%here method=1 is used, whereas all other calls of lyapunov_symm use method=2. The reason is that T and U are persistent, and input matrix A is the same, so using option 2 for all the rest of iterations spares a lot of computing time while not repeating Schur every time
indzeros = find(abs(Eyy0) < 1e-12); %find values that are numerical zero
- Eyy0(indzeros) = 0;
+ Eyy0(indzeros) = 0;
if useautocorr
sdy = sqrt(diag(Eyy0)); %theoretical standard deviation
+ sdy = sdy(stationary_vars);
sy = sdy*sdy'; %cross products of standard deviations
end
@@ -289,13 +296,13 @@ if ~no_identification_moments
else
dMOMENTS(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , jp) = dyn_vech(dEyy0); %focus only on VAROBS variables
end
- tmpEyyi = A*Ezz0*C' + B*Varinov*D';
+ tmpEyyi = A*Ezz0*C(stationary_vars,:)' + B*Varinov*D(stationary_vars,:)';
%we could distinguish between stderr and corr params, but this has no real speed effect as we multipliy with zeros
dtmpEyyi = dA(:,:,jp)*Ezz0*C' + A*dEzz0*C' + A*Ezz0*dC(:,:,jp)' + dB(:,:,jp)*Varinov*D' + B*dVarinov(:,:,jp)*D' + B*Varinov*dD(:,:,jp)';
Ai = eye(size(A,1)); %this is A^0
dAi = zeros(size(A,1),size(A,1)); %this is d(A^0)
for i = 1:nlags
- Eyyi = C*Ai*tmpEyyi;
+ Eyyi = C(stationary_vars,:)*Ai*tmpEyyi;
dEyyi = dC(:,:,jp)*Ai*tmpEyyi + C*dAi*tmpEyyi + C*Ai*dtmpEyyi;
if useautocorr
dEyyi = (dEyyi.*sy-dsy.*Eyyi)./(sy.*sy);
diff --git a/matlab/identification_numerical_objective.m b/matlab/identification_numerical_objective.m
index c7acf5e4b4d9960f20a6fc5e53ce886748a20eb1..b9af94a1ab83a92a135456a20fccf99ec7c4fef6 100644
--- a/matlab/identification_numerical_objective.m
+++ b/matlab/identification_numerical_objective.m
@@ -88,33 +88,53 @@ else
end
%% compute Kalman transition matrices and steady state with updated parameters
-[A, B, ~, ~, M, options, oo] = dynare_resolve(M, options, oo);
-ys = oo.dr.ys; %steady state of model variables in declaration order
-y0 = ys(oo.dr.order_var); %steady state of model variables in DR order
-
+[~,info,M,options,oo] = resol(0,M,options,oo);
+options = rmfield(options,'options_ident');
+oo.dr = pruned_state_space_system(M, options, oo.dr);
+A = oo.dr.pruned.A;
+B = oo.dr.pruned.B;
+C = oo.dr.pruned.C(indvar,:);
+D = oo.dr.pruned.D(indvar,:);
+Om_z = oo.dr.pruned.Om_z;
+Om_y = oo.dr.pruned.Om_y(indvar,indvar);
+Varinov = oo.dr.pruned.Varinov;
+obs_nbr = size(C,1);
%% out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is Iskrev (2010)'s J matrix.
if outputflag == 1
% Denote Ezz0 = E_t(z_t * z_t'), then the following Lyapunov equation defines the autocovariagrom: Ezz0 -A*Ezz*A' = B*Sig_e*B'
- Ezz0 = lyapunov_symm(A,B*M.Sigma_e*B',options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,[],options.debug);
- indzeros = find(abs(Ezz0) < 1e-12); %set small values to zero
- Ezz0(indzeros) = 0;
+ [Ezz0,u] = lyapunov_symm(A, Om_z, options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug);
+ stationary_vars = (1:size(C,1))';
+ if ~isempty(u)
+ x = abs(C*u);
+ stationary_vars = find(all(x < options.Schur_vec_tol,2));
+ end
+ Eyy0 = NaN*ones(obs_nbr,obs_nbr);
+ Eyy0(stationary_vars,stationary_vars) = C(stationary_vars,:)*Ezz0*C(stationary_vars,:)' + Om_y(stationary_vars,stationary_vars);
+ indzeros = find(abs(Eyy0) < 1e-12); %find values that are numerical zero
+ Eyy0(indzeros) = 0;
if useautocorr
- sy = sqrt(diag(Ezz0));
- sy = sy*sy';
+ sy = sqrt(diag(Ezz0)); %theoretical standard deviation
+ sy = sy(stationary_vars);
+ sy = sy*sy'; %cross products of standard deviations
sy0 = sy-diag(diag(sy))+eye(length(sy));
- Ezz0corr = Ezz0./sy0;
- out = dyn_vech(Ezz0corr(indvar,indvar)); %focus only on unique terms
+ Eyy0corr = NaN*ones(size(C,1),size(C,1));
+ Eyy0corr(stationary_vars,stationary_vars) = Eyy0./sy0;
+ out = dyn_vech(Eyy0corr); %focus only on unique terms
else
- out = dyn_vech(Ezz0(indvar,indvar)); %focus only on unique terms
+ out = dyn_vech(Eyy0); %focus only on unique terms
end
% compute autocovariances/autocorrelations of lagged observed variables
- for ii = 1:nlags
- Ezzii = A^(ii)*Ezz0;
+ tmpEyyi = A*Ezz0*C(stationary_vars,:)' + B*Varinov*D(stationary_vars,:)';
+ Ai = eye(size(A,1)); %this is A^0
+ for i = 1:nlags
+ Eyyi = NaN*ones(obs_nbr,obs_nbr);
+ Eyyi(stationary_vars,stationary_vars) = C(stationary_vars,:)*Ai*tmpEyyi;
if useautocorr
- Ezzii = Ezzii./sy;
+ Eyyi = Eyyi./sy;
end
- out = [out;vec(Ezzii(indvar,indvar))];
- end
+ out = [out;vec(Eyyi)];
+ Ai = Ai*A; %note that this is A^(i-1)
+ end
end
%% out = vec(g_omega). This is needed for Qu and Tkachenko (2012)'s G matrix.
@@ -122,15 +142,13 @@ if outputflag == 2
% This computes the spectral density g_omega where the interval [-pi;\pi] is discretized by grid_nbr points
freqs = (0 : pi/(grid_nbr/2):pi);% we focus only on positive values including the 0 frequency
tpos = exp( sqrt(-1)*freqs); %Fourier frequencies
- C = A(indvar,:);
- D = B(indvar,:);
IA = eye(size(A,1));
- var_nbr = length(indvar);
+ var_nbr = size(C,1);
out = zeros(var_nbr^2*length(freqs),1);
kk = 0;
for ig = 1:length(freqs)
Transferfct = D + C*((tpos(ig)*IA-A)\B);
- g_omega = (1/(2*pi))*(Transferfct*M.Sigma_e*Transferfct'); % note that ' is the conjugate transpose
+ g_omega = (1/(2*pi))*(Transferfct*Varinov*Transferfct'); % note that ' is the conjugate transpose
kk = kk+1;
out(1 + (kk-1)*var_nbr^2 : kk*var_nbr^2) = g_omega(:);
end
diff --git a/tests/identification/ident_unit_root/ident_unit_root_xfail.mod b/tests/identification/ident_unit_root/ident_unit_root_xfail.mod
index 45186c261b4e9ac266ba52f6c0af5685eb82c496..e8c2744de915ccc22168eca501cf433dc2b2eb8b 100644
--- a/tests/identification/ident_unit_root/ident_unit_root_xfail.mod
+++ b/tests/identification/ident_unit_root/ident_unit_root_xfail.mod
@@ -1,7 +1,7 @@
//Tests Identification command with ML and unit roots/diffuse filter option;
//Should not work because of observed unit root variable
-var y delta_y x z;
+var y x z delta_y;
varexo eps_x eps_z;