Fixed bug regarding non-stationary variables in pruned moments

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

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);

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    

tests/identification/ident_unit_root/ident_unit_root_xfail.mod

 //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;