diff --git a/matlab/+mom/check_irf_matching_file.m b/matlab/+mom/check_irf_matching_file.m
index 70eb5b37a662cebe94ed8f25c969caab77a28e51..3aa5db2726446db02de2cb6b6686b08a9806d208 100644
--- a/matlab/+mom/check_irf_matching_file.m
+++ b/matlab/+mom/check_irf_matching_file.m
@@ -6,7 +6,7 @@ function [irf_matching_file_name, irf_matching_file_path] = check_irf_matching_f
 % -------------------------------------------------------------------------
 % INPUTS
 % - irf_matching_file:  [string] user provided name (with possible path and extension)
-%                                of the MATLAB function that transforms model irfs
+%                                of the MATLAB function that transforms model IRFs
 % -------------------------------------------------------------------------
 % OUTPUTS
 % - irf_matching_file_name: [string] name of the MATLAB function (without extension)
diff --git a/matlab/+mom/default_option_mom_values.m b/matlab/+mom/default_option_mom_values.m
index 6504af5666c972f1c0a1c86e3584d9ff5ced4ffb..5ce9d34a94e491d76a27b9fb754a4d2c6dccf850 100644
--- a/matlab/+mom/default_option_mom_values.m
+++ b/matlab/+mom/default_option_mom_values.m
@@ -9,8 +9,8 @@ function options_mom_ = default_option_mom_values(options_mom_, options_, dname,
 % control of the estimation instead of possibly having to deal with options
 % chosen somewhere else in the mod file.
 % Note 2: we call a "mode" the minimum of the objective function, i.e.
-% the parameter vector that minimizes the distance between the moments/irfs
-% computed from the model and the moments/irfs computed from the data.
+% the parameter vector that minimizes the distance between the moments/IRFs
+% computed from the model and the moments/IRFs computed from the data.
 % -------------------------------------------------------------------------
 % INPUTS
 %  o options_mom_:           [structure]  all user-specified settings (from the method_of_moments command)
diff --git a/matlab/+mom/display_comparison_moments_irfs.m b/matlab/+mom/display_comparison_moments_irfs.m
index cdefd80c7f7c8520c7d4d228c377795c21a45322..da31a4e34a0fafb19826575b1f8a279062782a6a 100644
--- a/matlab/+mom/display_comparison_moments_irfs.m
+++ b/matlab/+mom/display_comparison_moments_irfs.m
@@ -1,7 +1,7 @@
 function display_comparison_moments_irfs(M_, options_mom_, data_moments, model_moments)
 % display_comparison_moments_irfs(M_, options_mom_, data_moments, model_moments)
 % -------------------------------------------------------------------------
-% Displays and saves to disk the comparison of the data moments/irfs and the model moments/irfs
+% Displays and saves to disk the comparison of the data moments/IRFs and the model moments/IRFs
 % -------------------------------------------------------------------------
 % INPUTS
 % M_:             [structure]  model information
@@ -40,7 +40,7 @@ function display_comparison_moments_irfs(M_, options_mom_, data_moments, model_m
 
 
 if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
-    titl = upper('Comparison of matched data irfs and model irfs');
+    titl = upper('Comparison of matched data IRFs and model IRFs');
     headers = {'IRF','Data','Model'};
     for jm = 1:size(M_.matched_irfs,1)
         labels{jm,1} = [M_.endo_names{M_.matched_irfs{jm,1}(1)} ' ' M_.exo_names{M_.matched_irfs{jm,1}(2)} ' (' num2str(M_.matched_irfs{jm,1}(3)) ')'];
diff --git a/matlab/+mom/graph_comparison_irfs.m b/matlab/+mom/graph_comparison_irfs.m
index 54aab693815cdc09c0db65df6aabcc208715f7bf..ee6d09e73bd09b1a1cc8a361cdff73b2b0f985ce 100644
--- a/matlab/+mom/graph_comparison_irfs.m
+++ b/matlab/+mom/graph_comparison_irfs.m
@@ -1,14 +1,14 @@
 function graph_comparison_irfs(matched_irfs,irf_model_varobs,varobs_id,irf_horizon,relative_irf,endo_names,exo_names,exo_names_tex,dname,fname,graph_format,TeX,nodisplay,figures_textwidth)
 % graph_comparison_irfs(matched_irfs,irf_model_varobs,varobs_id,irf_horizon,relative_irf,endo_names,exo_names,exo_names_tex,dname,fname,graph_format,TeX,nodisplay,figures_textwidth)
 % -------------------------------------------------------------------------
-% Plots and saves to disk the comparison of the selected data irfs and corresponding model irfs
+% Plots and saves to disk the comparison of the selected data IRFs and corresponding model IRfs
 % -------------------------------------------------------------------------
 % INPUTS
-% matched_irfs:      [matrix]   information on matched data irfs
-% irf_model_varobs:  [matrix]   model irfs for observable variables
+% matched_irfs:      [matrix]   information on matched data IRFs
+% irf_model_varobs:  [matrix]   model IRFs for observable variables
 % varobs_id:         [vector]   index for observable variables in endo_names
-% irf_horizon:       [scalar]   maximum horizon of irfs
-% relative_irf:      [boolean]  if true, plots normalized irfs
+% irf_horizon:       [scalar]   maximum horizon of IRFs
+% relative_irf:      [boolean]  if true, plots normalized IRFs
 % endo_names:        [cell]     names of endogenous variables
 % exo_names:         [cell]     names of exogenous variables
 % exo_names_tex:     [cell]     names of exogenous variables in latex
diff --git a/matlab/+mom/matched_irfs_blocks.m b/matlab/+mom/matched_irfs_blocks.m
index 0f91eed6ed7f92a03cef90722dd354dc4cb215c2..cce2d827816f5b17d1753072eb68aa0c223e078c 100644
--- a/matlab/+mom/matched_irfs_blocks.m
+++ b/matlab/+mom/matched_irfs_blocks.m
@@ -13,9 +13,9 @@ function [data_irfs, weight_mat, irf_index, max_irf_horizon] = matched_irfs_bloc
 % endo_names:          [cell array] list of endogenous variables
 % -------------------------------------------------------------------------
 % OUTPUT
-% data_irfs:           [matrix]     irfs for VAROBS as declared in matched_irfs block
-% weight_mat:          [matrix]     weighting matrix for irfs as declared in matched_irfs_weight block
-% irf_index:           [vector]     index for selecting specific irfs from full IRF matrix of observables
+% data_irfs:           [matrix]     IRFs for VAROBS as declared in matched_irfs block
+% weight_mat:          [matrix]     weighting matrix for IRFs as declared in matched_irfs_weight block
+% irf_index:           [vector]     index for selecting specific IRFs from full IRF matrix of observables
 % max_irf_horizon:     [scalar]     maximum IRF horizon as declared in matched_irfs block
 % -------------------------------------------------------------------------
 % This function is called by
@@ -42,8 +42,8 @@ function [data_irfs, weight_mat, irf_index, max_irf_horizon] = matched_irfs_bloc
 
 max_irf_horizon = max(cellfun(@(x) x(end), matched_irfs(:,1))); % get maximum IRF horizon
 % create full matrix where 1st dimension are IRF periods, 2nd dimension are variables as declared in VAROBS, 3rd dimension are shocks.
-data_irfs = nan(max_irf_horizon,obs_nbr,exo_nbr);
-% overwrite nan values if they are declared in matched_irfs block; remaining nan values will be later ignored in the matching
+data_irfs = NaN(max_irf_horizon,obs_nbr,exo_nbr);
+% overwrite NaN values if they are declared in matched_irfs block; remaining NaN values will be later ignored in the matching
 for jj = 1:size(matched_irfs,1)
     id_var       = matched_irfs{jj,1}(1);
     id_varobs    = find(varobs_id==id_var,1);
@@ -57,7 +57,7 @@ for jj = 1:size(matched_irfs,1)
     data_irfs(id_irf_period,id_varobs,id_shock) = irf_value;
 end
 % create (full) empirical weighting matrix
-weight_mat = eye(max_irf_horizon*obs_nbr*exo_nbr); % identity matrix by default: all irfs are equally important
+weight_mat = eye(max_irf_horizon*obs_nbr*exo_nbr); % identity matrix by default: all IRFs are equally important
 for jj = 1:size(matched_irfs_weight,1)
     id_var1 = matched_irfs_weight{jj,1}(1);  id_varobs1 = find(varobs_id==id_var1,1);  id_shock1 = matched_irfs_weight{jj,1}(2);  id_irf_period1 = matched_irfs_weight{jj,1}(3);
     id_var2 = matched_irfs_weight{jj,2}(1);  id_varobs2 = find(varobs_id==id_var2,1);  id_shock2 = matched_irfs_weight{jj,2}(2);  id_irf_period2 = matched_irfs_weight{jj,2}(3);
@@ -75,7 +75,7 @@ for jj = 1:size(matched_irfs_weight,1)
     weight_mat(idweight_mat1,idweight_mat2) = weight_mat_value;
     weight_mat(idweight_mat2,idweight_mat1) = weight_mat_value; % symmetry
 end
-% focus only on specified irfs
+% focus only on specified IRFs
 irf_index = find(~isnan(data_irfs));
 data_irfs = data_irfs(irf_index);
 weight_mat = weight_mat(irf_index,irf_index);
\ No newline at end of file
diff --git a/matlab/+mom/mode_compute_irf_matching.m b/matlab/+mom/mode_compute_irf_matching.m
index 4ae7de813c2c3ff49694967831e672baee85bae6..d7673b1261e98556933a4f29eff217a0e8555dbf 100644
--- a/matlab/+mom/mode_compute_irf_matching.m
+++ b/matlab/+mom/mode_compute_irf_matching.m
@@ -1,11 +1,11 @@
 function [xparam1, hessian_xparam1, fval, mom_verbose] = mode_compute_irf_matching(xparam0, hessian_xparam0, objective_function, doBayesianEstimation, weighting_info, data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
 % [xparam1, hessian_xparam1, fval, mom_verbose] = mode_compute_irf_matching(xparam0, hessian_xparam0, objective_function, doBayesianEstimation, weighting_info, data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
 % -------------------------------------------------------------------------
-% Computes the minimum of the objective function (distance between data irfs
-% and model irfs) for a sequence of optimizers.
+% Computes the minimum of the objective function (distance between data IRFs
+% and model IRFs) for a sequence of optimizers.
 % Note that we call a "mode" the minimum of the objective function, i.e.
-% the parameter vector that minimizes the distance between the irfs
-% computed from the model and the irfs computed from the data.
+% the parameter vector that minimizes the distance between the IRFs
+% computed from the model and the IRFs computed from the data.
 % -------------------------------------------------------------------------
 % INPUTS
 % xparam0:               [vector]       initialized parameters
diff --git a/matlab/+mom/objective_function.m b/matlab/+mom/objective_function.m
index bbef5e322f66934bc8d87cb4823532503bfa9216..fdcbd2365ee22d80c431346a020a90148f242e43 100644
--- a/matlab/+mom/objective_function.m
+++ b/matlab/+mom/objective_function.m
@@ -2,11 +2,11 @@ function [fval, info, exit_flag, df, junk_hessian, Q, model_moments, model_momen
 % [fval, info, exit_flag, df, junk_hessian, Q, model_moments, model_moments_params_derivs, irf_model_varobs] = objective_function(xparam, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
 % -------------------------------------------------------------------------
 % This function evaluates the objective function for method of moments estimation,
-% i.e. distance between model and data moments/irfs, possibly augmented with priors.
+% i.e. distance between model and data moments/IRFs, possibly augmented with priors.
 % -------------------------------------------------------------------------
 % INPUTS (same ones as in dsge_likelihood.m and dsge_var_likelihood.m)
 %  - xparam:               [vector]     current value of estimated parameters as returned by set_prior()
-%  - data_moments:         [vector]     data with moments/irfs to match (corresponds to dataset_ in likelihood-based estimation)
+%  - data_moments:         [vector]     data with moments/IRFs to match (corresponds to dataset_ in likelihood-based estimation)
 %  - weighting_info:       [structure]  storing information on weighting matrices (corresponds to dataset_info in likelihood-based estimation)
 %  - options_mom_:         [structure]  information about all settings (specified by the user, preprocessor, and taken from global options_)
 %  - M_                    [structure]  model information
@@ -27,7 +27,7 @@ function [fval, info, exit_flag, df, junk_hessian, Q, model_moments, model_momen
 %  - Q:                            [double]  value of the quadratic form of the moment difference
 %  - model_moments:                [vector]  model moments
 %  - model_moments_params_derivs:  [matrix]  analytical jacobian of the model moments wrt estimated parameters (currently for GMM only)
-%  - irf_model_varobs:             [matrix]  model irfs for observable variables (used for plotting matched irfs in mom.run)
+%  - irf_model_varobs:             [matrix]  model IRFs for observable variables (used for plotting matched IRfs in mom.run)
 % -------------------------------------------------------------------------
 % This function is called by
 %  o mom.run
@@ -78,12 +78,12 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
     if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
         if options_mom_.mom.vector_output == 1
             if options_mom_.mom.penalized_estimator
-                df = nan(options_mom_.mom.mom_nbr+length(xparam),length(xparam));
+                df = NaN(options_mom_.mom.mom_nbr+length(xparam),length(xparam));
             else
-                df = nan(options_mom_.mom.mom_nbr,length(xparam));
+                df = NaN(options_mom_.mom.mom_nbr,length(xparam));
             end
         else
-            df = nan(length(xparam),1);
+            df = NaN(length(xparam),1);
         end
     end
 end
@@ -226,7 +226,7 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
     scaled_shock_series(:,i_exo_var) = options_mom_.mom.shock_series(:,i_exo_var)*chol_S; % set non-zero entries
     % simulate series
     y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order);
-    % provide meaningful penalty if data is nan or inf
+    % provide meaningful penalty if data is NaN or Inf
     if any(any(isnan(y_sim))) || any(any(isinf(y_sim)))
         info(1)=180;
         info(4) = 0.1;
@@ -255,14 +255,14 @@ end
 
 
 %------------------------------------------------------------------------------
-% IRF_MATCHING using STOCH_SIMUL: Compute irfs given model solution and Cholesky
+% IRF_MATCHING using STOCH_SIMUL: Compute IRFs given model solution and Cholesky
 % decomposition on shock covariance matrix; this resembles the core codes in
 % stoch_simul
 %------------------------------------------------------------------------------
 if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL')
     cs = get_lower_cholesky_covariance(M_.Sigma_e,options_mom_.add_tiny_number_to_cholesky);
     irf_shocks_indx = getIrfShocksIndx(M_, options_mom_);
-    model_irf = nan(options_mom_.irf,M_.endo_nbr,M_.exo_nbr);
+    model_irf = NaN(options_mom_.irf,M_.endo_nbr,M_.exo_nbr);
     for i = irf_shocks_indx
         if options_mom_.order>1 && options_mom_.relative_irf % normalize shock to 0.01 before IRF generation for GIRFs; multiply with 100 later
             y_irf = irf(M_, options_mom_, dr, cs(:,i)./cs(i,i)/100, options_mom_.irf, options_mom_.drop, options_mom_.replic, options_mom_.order);
@@ -287,7 +287,7 @@ if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom
         end
         model_irf(:,:,i) = transpose(y_irf);
     end
-    % do transformations on model irfs if irf_matching_file is provided
+    % do transformations on model IRFs if irf_matching_file is provided
     if ~isempty(options_mom_.mom.irf_matching_file.name)
         [model_irf, check] = feval(str2func(options_mom_.mom.irf_matching_file.name), model_irf, M_, options_mom_, dr.ys);
         if check
diff --git a/matlab/+mom/print_info_on_estimation_settings.m b/matlab/+mom/print_info_on_estimation_settings.m
index 5d6ad6bfff8596d5cdc83acf3e10806e80742591..4085f1e0f1344f1b55aff9fe18ef7071b451d583 100644
--- a/matlab/+mom/print_info_on_estimation_settings.m
+++ b/matlab/+mom/print_info_on_estimation_settings.m
@@ -131,7 +131,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
     end
     fprintf('\n  - number of matched moments: %d', options_mom_.mom.mom_nbr);
 elseif strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
-    fprintf('\n  - number of matched irfs: %d', options_mom_.mom.mom_nbr);
+    fprintf('\n  - number of matched IRFs: %d', options_mom_.mom.mom_nbr);
 end
 fprintf('\n  - number of parameters: %d', number_of_estimated_parameters);
 fprintf('\n\n');
\ No newline at end of file
diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m
index fb1677931c5b17f69a86042f4ed4bea27d56030b..e970ccd278df04132573ec93e323377656b813d5 100644
--- a/matlab/+mom/run.m
+++ b/matlab/+mom/run.m
@@ -17,12 +17,12 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
 %  o Bayesian MCMC estimation
 %  o Display of results
 %    - GMM/SMM: J-Test and fit of moments
-%    - IRF_MATCHING: fit of irfs
+%    - IRF_MATCHING: fit of IRFs
 %  o Clean up
 % -------------------------------------------------------------------------
 % Note that we call a "mode" the minimum of the objective function, i.e.
-% the parameter vector that minimizes the distance between the moments/irfs
-% computed from the model and the moments/irfs computed from the data.
+% the parameter vector that minimizes the distance between the moments/IRFs
+% computed from the model and the moments/IRFs computed from the data.
 % -------------------------------------------------------------------------
 % This function is inspired by replication codes accompanied to the following papers:
 % GMM/SMM:
@@ -45,7 +45,7 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
 %                                                vars: matched_moments{:,1});
 %                                                lead/lags: matched_moments{:,2};
 %                                                powers: matched_moments{:,3};
-%                     o matched_irfs:          [cell] information about selected irfs to match in IRF_MATCHING estimation
+%                     o matched_irfs:          [cell] information about selected IRFs to match in IRF_MATCHING estimation
 %                     o matched_irfs_weights:  [cell] information about entries in weight matrix for an IRF_MATCHING estimation
 %  o options_mom_:   [structure] information about settings specified by the user
 % -------------------------------------------------------------------------
@@ -238,16 +238,16 @@ if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
         error('method_of_moments: Something wrong while computing inv(W), check your weighting matrix!');
     end
     if any(isnan(oo_.mom.weighting_info.Winv(:))) || any(isinf(oo_.mom.weighting_info.Winv(:)))
-        error('method_of_moments: There are nan or inf in inv(W), check your weighting matrix!');
+        error('method_of_moments: There are NaN or Inf values in inv(W), check your weighting matrix!');
     end
-    % compute log determinant of inverse of weighting matrix in a robust way to avoid inf/nan
+    % compute log determinant of inverse of weighting matrix in a robust way to avoid Inf or NaN
     try
         oo_.mom.weighting_info.Winv_logdet = 2*sum(log(diag(chol(oo_.mom.weighting_info.Winv))));
     catch
         error('method_of_moments: Something wrong while computing log(det(inv(W))), check your weighting matrix!');
     end
     if any(isnan(oo_.mom.weighting_info.Winv_logdet(:))) || any(isinf(oo_.mom.weighting_info.Winv_logdet(:)))
-        error('method_of_moments: There are nan or inf in log(det(inv(W))), check your weighting matrix!');
+        error('method_of_moments: There are NaN or Inf values in log(det(inv(W))), check your weighting matrix!');
     end
     options_mom_.mom.mom_nbr = length(options_mom_.mom.irfIndex);
 end
@@ -606,12 +606,12 @@ else
         end
     else
         variances = bayestopt_.p2.*bayestopt_.p2;
-        id_inf = isinf(variances);
-        variances(id_inf) = 1;
+        id_Inf = isinf(variances);
+        variances(id_Inf) = 1;
         invhess = options_mom_.mh_posterior_mode_estimation*diag(variances);
         xparam1 = bayestopt_.p5;
-        id_nan = isnan(xparam1);
-        xparam1(id_nan) = bayestopt_.p1(id_nan);
+        id_NaN = isnan(xparam1);
+        xparam1(id_NaN) = bayestopt_.p1(id_NaN);
         outside_bound_pars=find(xparam1 < BoundsInfo.lb | xparam1 > BoundsInfo.ub);
         xparam1(outside_bound_pars) = bayestopt_.p1(outside_bound_pars);
     end
@@ -802,7 +802,7 @@ elseif strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
         mom.graph_comparison_irfs(M_.matched_irfs,oo_.mom.irf_model_varobs,options_mom_.varobs_id,options_mom_.irf,options_mom_.relative_irf,M_.endo_names,M_.exo_names,M_.exo_names_tex,M_.dname,M_.fname,options_mom_.graph_format,options_mom_.TeX,options_mom_.nodisplay,options_mom_.figures.textwidth)
     end
 end
-% display comparison of model moments/irfs and data moments/irfs
+% display comparison of model moments/IRFs and data moments/IRFs
 mom.display_comparison_moments_irfs(M_, options_mom_, oo_.mom.data_moments, oo_.mom.model_moments);
 % save results to _mom_results.mat
 save([M_.dname filesep 'method_of_moments' filesep M_.fname '_mom_results.mat'], 'oo_', 'options_mom_', 'M_', 'estim_params_', 'bayestopt_');
diff --git a/matlab/+mom/standard_errors.m b/matlab/+mom/standard_errors.m
index 55c64d0d28c3c2b2c22551f533202e5fcd0ba06e..c8336d68faf416b9eaeab368365a12ed8affc654 100644
--- a/matlab/+mom/standard_errors.m
+++ b/matlab/+mom/standard_errors.m
@@ -12,7 +12,7 @@ function [stderr_values, asympt_cov_mat] = standard_errors(xparam, objective_fun
 %  - model_moments:        [vector]     model moments
 %  - model_moments_params_derivs:  [matrix]  analytical jacobian of the model moments wrt estimated parameters (currently for GMM only)
 %  - m_data                [matrix]     selected empirical moments at each point in time
-%  - data_moments:         [vector]     data with moments/irfs to match
+%  - data_moments:         [vector]     data with moments/IRFs to match
 %  - weighting_info:       [structure]  storing information on weighting matrices
 %  - options_mom_:         [structure]  information about all settings (specified by the user, preprocessor, and taken from global options_)
 %  - M_                    [structure]  model information