Merge branch 'MoM_fixes' into 'master'

Method of Moments: Add GMM analytical standard errors and some other minor fixes

See merge request Dynare/dynare!1791

matlab/method_of_moments/method_of_moments.m

@@ -86,21 +86,12 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
 % =========================================================================
 
 %% TO DO LIST
-% - [ ] why does lsqnonlin take less time in Andreasen toolbox?
-% - [ ] test user-specified weightning matrix
-% - [ ] which qz_criterium value?
-% - [ ] document that in method_of_moments_data_moments.m NaN are replaced by mean of moment
 % - [ ] add IRF matching
-% - [ ] test estimated_params_bounds block
-% - [ ] test what happens if all parameters will be estimated but some/all are not calibrated
-% - [ ] speed up lyapunov equation by using doubling with old initial values
-% - [ ] check smm at order > 3 without pruning
-% - [ ] provide option to use analytical derivatives to compute std errors (similar to what we already do in identification)
-% - [ ] add Bayesian GMM/SMM estimation
-% - [ ] useautocorr
-% - [ ] do we need dirname?
-% - [ ] decide on default weighting matrix scheme, I would propose 2 stage with Diagonal of optimal matrix
-% - [ ] check smm with product moments greater than 2
+% - [ ] speed up pruned_state_space_system (by using doubling with old initial values, hardcoding zeros, other "tricks" used in e.g. nlma)
+% - [ ] add option to use autocorrelations (we have useautocorr in identification toolbox already)
+% - [ ] SMM with extended path
+% - [ ] deal with measurement errors (once @wmutschl has implemented this in identification toolbox)
+% - [ ] improve check for duplicate moments by using the cellfun and unique functions
 % -------------------------------------------------------------------------
 % Step 0: Check if required structures and options exist
 % -------------------------------------------------------------------------
@@ -125,7 +116,7 @@ else
     options_mom_.loglinear = false;
 end
 
-fprintf('\n==== Method of Moments (%s) Estimation ====\n\n',options_mom_.mom.mom_method)
+fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom_method)
 
 % -------------------------------------------------------------------------
 % Step 1a: Prepare options_mom_ structure
@@ -147,10 +138,10 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
     options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix_scaling_factor',1);    % scaling of weighting matrix
     options_mom_.mom = set_default_option(options_mom_.mom,'se_tolx',1e-5);                         % step size for numerical computation of standard errors
     options_mom_ = set_default_option(options_mom_,'order',1);                                      % order of Taylor approximation in perturbation
-    options_mom_ = set_default_option(options_mom_,'pruning',true);                                 % use pruned state space system at higher-order
+    options_mom_ = set_default_option(options_mom_,'pruning',false);                                % use pruned state space system at higher-order
     % Checks for perturbation order
     if options_mom_.order < 1
-        error('method_of_moments:: The order of the Taylor approximation cannot be 0!')
+        error('method_of_moments: The order of the Taylor approximation cannot be 0!')
     end
 end
 if strcmp(options_mom_.mom.mom_method,'SMM')
@@ -169,11 +160,15 @@ 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
+    if options_mom_.order > 3
+        error('method_of_moments: perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM.\n');
+    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
-options_mom_ = set_default_option(options_mom_,'dirname',M_.dname);    % directory in which to store estimation output
 options_mom_ = set_default_option(options_mom_,'graph_format','eps');  % specify the file format(s) for graphs saved to disk
 options_mom_ = set_default_option(options_mom_,'nodisplay',false);     % do not display the graphs, but still save them to disk
 options_mom_ = set_default_option(options_mom_,'nograph',false);       % do not create graphs (which implies that they are not saved to the disk nor displayed)
@@ -226,8 +221,11 @@ options_mom_ = set_default_option(options_mom_,'lyapunov_fixed_point_tol',1e-10)
 options_mom_ = set_default_option(options_mom_,'lyapunov_doubling_tol',1e-16);          % convergence criterion used in the doubling algorithm
 options_mom_ = set_default_option(options_mom_,'sylvester_fp',false);                   % determines whether to use fixed point algorihtm to solve Sylvester equation (gensylv_fp), faster for large scale models
 options_mom_ = set_default_option(options_mom_,'sylvester_fixed_point_tol',1e-12);      % convergence criterion used in the fixed point Sylvester solver
-options_mom_ = set_default_option(options_mom_,'qz_criterium',1-1e-6);                  % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [IS THIS CORRET @wmutschl]
+options_mom_ = set_default_option(options_mom_,'qz_criterium',1-1e-6);                  % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
+                                                                                        % if there are no unit roots one can use 1.0 (or slightly below) which we set as default; if they are possible, you may have have multiple unit roots and the accuracy decreases when computing the eigenvalues in lyapunov_symm
+                                                                                        % Note that unit roots are only possible at first-order, at higher order we set it to 1 in pruned_state_space_system and focus only on stationary observables.
 options_mom_ = set_default_option(options_mom_,'qz_zero_threshold',1e-6);               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
+options_mom_ = set_default_option(options_mom_,'schur_vec_tol',1e-11);                  % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix.
 if options_mom_.order > 2
     fprintf('Dynare will use ''k_order_solver'' as the order>2\n');
     options_mom_.k_order_solver = true;
@@ -326,6 +324,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
 
@@ -362,8 +361,6 @@ end
 % -------------------------------------------------------------------------
 % Step 2: Checks and transformations for matched moments structure (preliminary)
 % -------------------------------------------------------------------------
-% Note that we do not have a preprocessor interface yet for this, so this
-% will need much improvement later on. @wmutschl
 
 % Initialize indices
 options_mom_.mom.index.E_y       = false(options_mom_.obs_nbr,1);                      %unconditional first order product moments
@@ -423,8 +420,6 @@ for jm=1:size(M_.matched_moments,1)
     end
 end
 
-
-% @wmutschl: add check for duplicate moments by using the cellfun and unique functions
 %Remove duplicate elements
 UniqueMomIdx = [nonzeros(options_mom_.mom.index.E_y_pos); nonzeros(tril(options_mom_.mom.index.E_yy_pos)); nonzeros(options_mom_.mom.index.E_yyt_pos)];
 DuplicateMoms = setdiff(1:size(M_.matched_moments,1),UniqueMomIdx);
@@ -614,16 +609,14 @@ if ~isempty(dataset_)
     options_mom_.nobs = dataset_.nobs;
 end
 
-% provide info on missing observations
-if any(any(isnan(dataset_.data)))
-    fprintf('missing observations will be replaced by the sample mean of the corresponding moment')
-end
-
 % Check length of data for estimation of second moments
 if options_mom_.ar > options_mom_.nobs+1
     error('method_of_moments: Data set is too short to compute second moments');
 end
 
+% Provide info on data moments handling
+fprintf('Computing data moments. Note that NaN values in the moments (due to leads and lags or missing data) are replaced by the mean of the corresponding moment\n');
+
 % Get data moments for the method of moments
 [oo_.mom.data_moments, oo_.mom.m_data] = method_of_moments_data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_);
 
@@ -651,6 +644,9 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
     end
     options_mom_.mom.shock_series = temp_shocks;
     options_mom_.mom.ME_shock_series = temp_shocks_ME;
+    if options_mom_.k_order_solver && ~options_mom_.pruning % dynare++ routines will be called in simult_.m, store some additional stuff
+        options_mom_.DynareRandomStreams.seed = options_mom_.mom.seed;
+    end
 end
 
 % -------------------------------------------------------------------------
@@ -688,6 +684,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 +786,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 +868,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

matlab/method_of_moments/method_of_moments_check_plot.m

@@ -39,11 +39,11 @@ end
 
 [nbplt,nr,nc,lr,lc,nstar] = pltorg(length(xparam));
 
-if ~exist([M_.fname filesep 'graphs'],'dir')
-    mkdir(M_.fname,'graphs');
+if ~exist([M_.dname filesep 'graphs'],'dir')
+    mkdir(M_.dname,'graphs');
 end
 if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
-    fidTeX = fopen([M_.fname, '/graphs/', M_.fname '_MoMCheckPlots.tex'],'w');
+    fidTeX = fopen([M_.dname, '/graphs/', M_.fname '_MoMCheckPlots.tex'],'w');
     fprintf(fidTeX,'%% TeX eps-loader file generated by method_of_moments_check_plot.m (Dynare).\n');
     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
     fprintf(fidTeX,' \n');
@@ -165,12 +165,12 @@ for plt = 1:nbplt
         text(0.25,0.5,'log-post')
         text(0.69,0.5,'log-lik kernel')
     end
-    dyn_saveas(hh,[M_.fname, '/graphs/', M_.fname '_MoMCheckPlots' int2str(plt) ],options_.nodisplay,options_.graph_format);
+    dyn_saveas(hh,[M_.dname, '/graphs/', M_.fname '_MoMCheckPlots' int2str(plt) ],options_.nodisplay,options_.graph_format);
     if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
         % TeX eps loader file
         fprintf(fidTeX,'\\begin{figure}[H]\n');
         fprintf(fidTeX,'\\centering \n');
-        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_MoMCheckPlots%s}\n',options_.figures.textwidth*min(k/nc,1),[M_.fname, '/graphs/',M_.fname],int2str(plt));
+        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_MoMCheckPlots%s}\n',options_.figures.textwidth*min(k/nc,1),[M_.dname, '/graphs/',M_.fname],int2str(plt));
         fprintf(fidTeX,'\\caption{Method of Moments check plots.}');
         fprintf(fidTeX,'\\label{Fig:MoMCheckPlots:%s}\n',int2str(plt));
         fprintf(fidTeX,'\\end{figure}\n');
@@ -181,5 +181,4 @@ if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
     fclose(fidTeX);
 end
 
-OutputDirectoryName = CheckPath('check_plot',M_.dname);
-save([OutputDirectoryName '/MoM_check_plot_data.mat'],'mcheck');
+save([M_.dname filesep 'graphs' filesep M_.fname '_MoMCheckPlots_data.mat'],'mcheck');

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

matlab/method_of_moments/method_of_moments_standard_errors.m

@@ -57,30 +57,46 @@ 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_);
-
-    % 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)
+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
+    idx_nan = find(any(isnan(D)));
+    if any(idx_nan)
+        for i = idx_nan            
+             fprintf('No standard errors available for parameter %s\n',get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_))
+        end        
+        warning('There are NaN in the analytical Jacobian of Moments. Check your bounds and/or priors, or use a different optimizer.')
         Asympt_Var = NaN(length(xparam),length(xparam));
         SE_values = NaN(length(xparam),1);
         return
     end
+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_);
+
+        % 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
 
 T = options_mom_.nobs; %Number of observations

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