diff --git a/matlab/method_of_moments/method_of_moments.m b/matlab/method_of_moments/method_of_moments.m
index c19a50ccf..6417ef22d 100644
--- a/matlab/method_of_moments/method_of_moments.m
+++ b/matlab/method_of_moments/method_of_moments.m
@@ -169,7 +169,9 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
         fprintf('GMM at higher order only works with pruning, so we set pruning option to 1.\n');
         options_mom_.pruning = true;
     end
+    options_mom_.mom = set_default_option(options_mom_.mom,'analytic_standard_errors',false);       % compute standard errors numerically (0) or analytically (1). Analytical derivatives are only available for GMM.
 end
+options_mom_.mom.compute_derivs = false;% flag to compute derivs in objective function (needed for analytic standard errors with GMM)
 
     
 % General options that can be set by the user in the mod file, otherwise default values are provided
@@ -326,6 +328,7 @@ options_mom_.solveopt         = options_.solveopt;
 options_mom_.gradient_method  = options_.gradient_method;
 options_mom_.gradient_epsilon = options_.gradient_epsilon;
 options_mom_.analytic_derivation = 0;
+options_mom_.analytic_derivation_mode = 0; % needed by get_perturbation_params_derivs.m, ie use efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012)
 
 options_mom_.vector_output= false;           % specifies whether the objective function returns a vector
 
@@ -688,6 +691,11 @@ if isfield(estim_params_,'param_vals') && ~isempty(estim_params_.param_vals)
         fprintf('This will override parameter values and may lead to wrong results.\n')
         fprintf('Check whether this is really intended.\n')
         warning('The steady state file internally changes the values of the estimated parameters.')
+        if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
+            fprintf('For analytical standard errors, the parameter-Jacobians of the dynamic model and of the steady-state will be computed numerically\n'),
+            fprintf('(re-set options_mom_.analytic_derivation_mode= -2)'),
+            options_mom_.analytic_derivation_mode= -2;
+        end
     end
 end
 
@@ -785,8 +793,14 @@ fprintf('\n  - perturbation order:        %d', options_mom_.order)
 if options_mom_.order > 1 && options_mom_.pruning
     fprintf(' (with pruning)')
 end
+if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
+    fprintf('\n  - standard errors:           analytic derivatives');
+else
+    fprintf('\n  - standard errors:           numerical derivatives');
+end
 fprintf('\n  - number of matched moments: %d', options_mom_.mom.mom_nbr);
-fprintf('\n  - number of parameters:      %d\n\n', length(xparam0));
+fprintf('\n  - number of parameters:      %d', length(xparam0));
+fprintf('\n\n');
 
 % -------------------------------------------------------------------------
 % Step 7b: Iterated method of moments estimation
@@ -861,8 +875,12 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
     options_mom_.vector_output = false;    
     % Update M_ and DynareResults (in particular to get oo_.mom.model_moments)    
     M_ = set_all_parameters(xparam1,estim_params_,M_);
+    if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
+        options_mom_.mom.compute_derivs = true; % for GMM we compute derivatives analytically in the objective function with this flag        
+    end
     [fval, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
-    % Compute Standard errors
+    options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization
+    
     SE = method_of_moments_standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Woptflag);
     
     % Store results in output structure
diff --git a/matlab/method_of_moments/method_of_moments_objective_function.m b/matlab/method_of_moments/method_of_moments_objective_function.m
index 8a15ed749..b97c7725f 100644
--- a/matlab/method_of_moments/method_of_moments_objective_function.m
+++ b/matlab/method_of_moments/method_of_moments_objective_function.m
@@ -109,7 +109,32 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
     %--------------------------------------------------------------------------
     % 3. Set up pruned state-space system and compute model moments
     %--------------------------------------------------------------------------
-    pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.dr.obs_var, options_mom_.ar, 0, 0);
+    if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors        
+        indpmodel = []; %initialize index for model parameters
+        if ~isempty(estim_params_.param_vals)
+            indpmodel = estim_params_.param_vals(:,1); %values correspond to parameters declaration order, row number corresponds to order in estimated_params
+        end
+        indpstderr=[]; %initialize index for stderr parameters
+        if ~isempty(estim_params_.var_exo)
+            indpstderr = estim_params_.var_exo(:,1); %values correspond to varexo declaration order, row number corresponds to order in estimated_params
+        end
+        indpcorr=[]; %initialize matrix for corr paramters
+        if ~isempty(estim_params_.corrx)
+            indpcorr = estim_params_.corrx(:,1:2); %values correspond to varexo declaration order, row number corresponds to order in estimated_params
+        end
+        if estim_params_.nvn || estim_params_.ncn %nvn is number of stderr parameters and ncn is number of corr parameters of measurement innovations as declared in estimated_params
+            error('Analytic computation of standard errrors does not (yet) support measurement errors.\nInstead, define them explicitly as varexo and provide measurement equations in the model definition.\nAlternatively, use numerical standard errors.')
+        end
+        modparam_nbr = estim_params_.np;        % number of model parameters as declared in estimated_params
+        stderrparam_nbr = estim_params_.nvx;    % number of stderr parameters
+        corrparam_nbr = estim_params_.ncx;      % number of corr parameters
+        totparam_nbr = stderrparam_nbr+corrparam_nbr+modparam_nbr;
+        dr.derivs = get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices
+        oo_.mom.model_moments_params_derivs = NaN(options_mom_.mom.mom_nbr,totparam_nbr);
+        pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.dr.obs_var, options_mom_.ar, 0, 1);
+    else
+        pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.dr.obs_var, options_mom_.ar, 0, 0);
+    end
     
     oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
     offset = 0;
@@ -118,6 +143,9 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
         E_y = pruned_state_space.E_y;
         E_y_nbr = nnz(options_mom_.mom.index.E_y);
         oo_.mom.model_moments(offset+1:E_y_nbr,1) = E_y(options_mom_.mom.index.E_y);
+        if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
+            oo_.mom.model_moments_params_derivs(offset+1:E_y_nbr,:) = pruned_state_space.dE_y(options_mom_.mom.index.E_y,:);
+        end
         offset = offset + E_y_nbr;
     end
     % Second moments
@@ -125,22 +153,47 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
     if isfield(options_mom_.mom.index,'E_yy') && nnz(options_mom_.mom.index.E_yy) > 0
         if options_mom_.prefilter
             E_yy = pruned_state_space.Var_y;
+            if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
+                dE_yy = pruned_state_space.dVar_y;
+            end            
         else
             E_yy = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y';
+            if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
+                dE_yy = pruned_state_space.dVar_y;
+                for jp=1:totparam_nbr
+                    dE_yy(:,:,jp) = dE_yy(:,:,jp) + pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y' + pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)';
+                end
+            end
         end
         E_yy_nbr = nnz(tril(options_mom_.mom.index.E_yy));
         oo_.mom.model_moments(offset+(1:E_yy_nbr),1) = E_yy(tril(options_mom_.mom.index.E_yy));
+        if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
+            oo_.mom.model_moments_params_derivs(offset+(1:E_yy_nbr),:) = reshape(dE_yy(repmat(tril(options_mom_.mom.index.E_yy),[1 1 totparam_nbr])),E_yy_nbr,totparam_nbr);
+        end
         offset = offset + E_yy_nbr;
     end
     % Lead/lags covariance
     if isfield(options_mom_.mom.index,'E_yyt') && nnz(options_mom_.mom.index.E_yyt) > 0
         if options_mom_.prefilter
             E_yyt = pruned_state_space.Var_yi;
+            if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
+                dE_yyt = pruned_state_space.dVar_yi;
+            end
         else
             E_yyt = pruned_state_space.Var_yi + repmat(pruned_state_space.E_y*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)]);
+            if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
+                dE_yyt = pruned_state_space.dVar_yi;
+                for jp=1:totparam_nbr
+                    dE_yyt(:,:,:,jp) = dE_yyt(:,:,:,jp) + repmat(pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)])...
+                                                        + repmat(pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)',[1 1 size(pruned_state_space.Var_yi,3)]);
+                end
+            end
         end
         E_yyt_nbr = nnz(options_mom_.mom.index.E_yyt);
         oo_.mom.model_moments(offset+(1:E_yyt_nbr),1) = E_yyt(options_mom_.mom.index.E_yyt);
+        if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
+            oo_.mom.model_moments_params_derivs(offset+(1:E_yyt_nbr),:) = reshape(dE_yyt(repmat(options_mom_.mom.index.E_yyt,[1 1 1 totparam_nbr])),E_yyt_nbr,totparam_nbr);
+        end
     end
 
 elseif strcmp(options_mom_.mom.mom_method,'SMM')
diff --git a/matlab/method_of_moments/method_of_moments_standard_errors.m b/matlab/method_of_moments/method_of_moments_standard_errors.m
index 888ccf01f..5084d47d0 100644
--- a/matlab/method_of_moments/method_of_moments_standard_errors.m
+++ b/matlab/method_of_moments/method_of_moments_standard_errors.m
@@ -57,29 +57,35 @@ dim_params   = size(xparam,1);
 D            = zeros(num_mom,dim_params);
 eps_value    = options_mom_.mom.se_tolx;
 
-for i=1:dim_params
-    %Positive step
-    xparam_eps_p      = xparam;
-    xparam_eps_p(i,1) = xparam_eps_p(i) + eps_value;
-    [~, info_p, ~, ~,~, oo__p] = feval(objective_function, xparam_eps_p, Bounds, oo_, estim_params_, M_, options_mom_);
-    
-    % Negative step
-    xparam_eps_m      = xparam;
-    xparam_eps_m(i,1) = xparam_eps_m(i) - eps_value;
-    [~, info_m,  ~, ~,~, oo__m] = feval(objective_function, xparam_eps_m, Bounds, oo_, estim_params_, M_, options_mom_);
+if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
+    fprintf('\nComputing standard errors using analytical derivatives of moments\n');
+    D = oo_.mom.model_moments_params_derivs; %already computed in objective function via get_perturbation_params.m
+else    
+    fprintf('\nComputing standard errors using numerical derivatives of moments\n');
+    for i=1:dim_params
+        %Positive step
+        xparam_eps_p      = xparam;
+        xparam_eps_p(i,1) = xparam_eps_p(i) + eps_value;
+        [~, info_p, ~, ~,~, oo__p] = feval(objective_function, xparam_eps_p, Bounds, oo_, estim_params_, M_, options_mom_);
 
-    % The Jacobian:
-    if nnz(info_p)==0 && nnz(info_m)==0
-        D(:,i) = (oo__p.mom.model_moments - oo__m.mom.model_moments)/(2*eps_value);
-    else
-        problpar = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_);
-        message_p = get_error_message(info_p, options_mom_);
-        message_m = get_error_message(info_m, options_mom_);        
-        
-        warning('method_of_moments:info','Cannot compute the Jacobian for parameter %s - no standard errors available\n %s %s\nCheck your bounds and/or priors, or use a different optimizer.\n',problpar, message_p, message_m)
-        Asympt_Var = NaN(length(xparam),length(xparam));
-        SE_values = NaN(length(xparam),1);
-        return
+        % Negative step
+        xparam_eps_m      = xparam;
+        xparam_eps_m(i,1) = xparam_eps_m(i) - eps_value;
+        [~, info_m,  ~, ~,~, oo__m] = feval(objective_function, xparam_eps_m, Bounds, oo_, estim_params_, M_, options_mom_);
+
+        % The Jacobian:
+        if nnz(info_p)==0 && nnz(info_m)==0
+            D(:,i) = (oo__p.mom.model_moments - oo__m.mom.model_moments)/(2*eps_value);
+        else
+            problpar = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_);
+            message_p = get_error_message(info_p, options_mom_);
+            message_m = get_error_message(info_m, options_mom_);        
+
+            warning('method_of_moments:info','Cannot compute the Jacobian for parameter %s - no standard errors available\n %s %s\nCheck your bounds and/or priors, or use a different optimizer.\n',problpar, message_p, message_m)
+            Asympt_Var = NaN(length(xparam),length(xparam));
+            SE_values = NaN(length(xparam),1);
+            return
+        end
     end
 end
 
diff --git a/tests/estimation/method_of_moments/AnScho_MoM.mod b/tests/estimation/method_of_moments/AnScho_MoM.mod
index 2f9c2c24d..83cec9136 100644
--- a/tests/estimation/method_of_moments/AnScho_MoM.mod
+++ b/tests/estimation/method_of_moments/AnScho_MoM.mod
@@ -25,7 +25,7 @@
 @#define estimParams = 1
 
 % Note that we set the numerical optimization tolerance levels very large to speed up the testsuite
-@#define optimizer = 4
+@#define optimizer = 13
 
 var c p R g y z INFL INT YGR;
 varexo e_r e_g e_z;
@@ -248,7 +248,7 @@ end
         %, optim = ('TolFun', 1e-5
         %           ,'TolX', 1e-6
         %          )    % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
-        %, silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
+        , silent_optimizer                  % run minimization of moments distance silently without displaying results or saving files in between
         % , tolf = 1e-5                       % convergence criterion on function value for numerical differentiation
         % , tolx = 1e-6                       % convergence criterion on funciton input for numerical differentiation
         % 
@@ -267,6 +267,9 @@ end
         % , sylvester_fixed_point_tol = 1e-12      % convergence criterion used in the fixed point Sylvester solver
         % , qz_criterium = 0.999999                % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [IS THIS CORRET @wmutschl]
         % , qz_zero_threshold = 1e-6               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
+        @#if mommethod == "GMM"
+        , analytic_standard_errors
+        @#endif
     );
 @#endfor