Fix wrong computations when identification(useautocorr)

(cherry picked from commit 347be1f8)

matlab/get_identification_jacobians.m

@@ -241,12 +241,11 @@ if order == 1
 end
 
 %% Compute dMOMENTS
-if ~no_identification_moments    
-    if useautocorr
-        E_yy  = pruned.Corr_y;  dE_yy  = pruned.dCorr_y;
+if ~no_identification_moments
+    E_yy  = pruned.Var_y;  dE_yy  = pruned.dVar_y;
+    if useautocorr        
         E_yyi = pruned.Corr_yi; dE_yyi = pruned.dCorr_yi;
-    else
-        E_yy  = pruned.Var_y;   dE_yy  = pruned.dVar_y;
+    else        
         E_yyi = pruned.Var_yi;  dE_yyi = pruned.dVar_yi;
     end
     MOMENTS = [MEAN; dyn_vech(E_yy)];

matlab/identification_numerical_objective.m

@@ -82,11 +82,7 @@ pruned = pruned_state_space_system(M, options, oo.dr, indvar, nlags, useautocorr
 
 %% 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    
-    if useautocorr
-        out = dyn_vech(pruned.Corr_y);
-    else
-        out = dyn_vech(pruned.Var_y);
-    end    
+    out = dyn_vech(pruned.Var_y);
     for i = 1:nlags
         if useautocorr
             out = [out;vec(pruned.Corr_yi(:,:,i))];

matlab/pruned_state_space_system.m

@@ -1032,9 +1032,9 @@ end
 indzeros = find(abs(Var_y) < 1e-12); %find values that are numerical zero
 Var_y(indzeros) = 0;
 if useautocorr
-    sy = sqrt(diag(Var_y)); %theoretical standard deviation
-    sy = sy(stationary_vars);
-    sy = sy*sy';            %cross products of standard deviations
+    sdy = sqrt(diag(Var_y)); %theoretical standard deviation
+    sdy = sdy(stationary_vars);
+    sy = sdy*sdy';           %cross products of standard deviations
     Corr_y = NaN*ones(y_nbr,y_nbr);
     Corr_y(stationary_vars,stationary_vars) = Var_y(stationary_vars,stationary_vars)./sy;
     Corr_yi = NaN*ones(y_nbr,y_nbr,nlags);
@@ -1060,10 +1060,9 @@ if compute_derivs
         dVar_y_tmp(indzeros) = 0;        
         dVar_y(stationary_vars,stationary_vars,jpV) = dVar_y_tmp;
         if useautocorr
-            %is this correct?[@wmutschl]
-            dsy = 1/2./sy.*diag(dVar_y(:,:,jpV));
+            dsy = 1/2./sdy.*diag(dVar_y(:,:,jpV));
             dsy = dsy(stationary_vars);
-            dsy = dsy*sy'+sy*dsy';
+            dsy = dsy*sdy'+sdy*dsy';
             dCorr_y(stationary_vars,stationary_vars,jpV) = (dVar_y(stationary_vars,stationary_vars,jpV).*sy-dsy.*Var_y(stationary_vars,stationary_vars))./(sy.*sy);
             dCorr_y(stationary_vars,stationary_vars,jpV) = dCorr_y(stationary_vars,stationary_vars,jpV)-diag(diag(dCorr_y(stationary_vars,stationary_vars,jpV)))+diag(diag(dVar_y(stationary_vars,stationary_vars,jpV)));
         end
@@ -1138,6 +1137,9 @@ if compute_derivs
                                                                 + dD(stationary_vars,:,jpVi)*E_inovzlagi*C(stationary_vars,:)' + D(stationary_vars,:)*dE_inovzlagi_jpVi*C(stationary_vars,:)' + D(stationary_vars,:)*E_inovzlagi*dC(stationary_vars,:,jpVi)';
             end
             if useautocorr
+                dsy = 1/2./sdy.*diag(dVar_y(:,:,jpVi));
+                dsy = dsy(stationary_vars);
+                dsy = dsy*sdy'+sdy*dsy';
                 dCorr_yi(stationary_vars,stationary_vars,i,jpVi) = (dVar_yi(stationary_vars,stationary_vars,i,jpVi).*sy-dsy.*Var_yi(stationary_vars,stationary_vars,i))./(sy.*sy);                
             end
             dAi_jpVi = dAi_jpVi*A + Ai*dA(:,:,jpVi);