diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/UnivariateSpectralDensity.m
index 5895c02e09804cd1a3d1becb9722af0810c59264..b46bdfd9c6566063af7e13014ada2bc5850dd2f2 100644
--- a/matlab/UnivariateSpectralDensity.m
+++ b/matlab/UnivariateSpectralDensity.m
@@ -150,7 +150,7 @@ else
     warning on MATLAB:dividebyzero
 end
 
-if options_.nograph == 0
+if ~options_.nograph
     if ~exist(M_.fname, 'dir')
         mkdir('.',M_.fname);
     end
diff --git a/matlab/check.m b/matlab/check.m
index 219ce1868847f6c5a96da8e8507182b8ada7ab6e..fba5bb8cd4f970e985957d57f8c5efae0f7cc897 100644
--- a/matlab/check.m
+++ b/matlab/check.m
@@ -66,7 +66,7 @@ if (nyf == dr.edim) && (dr.full_rank)
     result = 1;
 end
 
-if options.noprint == 0
+if ~options.noprint
     skipline()
     disp('EIGENVALUES:')
     disp(sprintf('%16s %16s %16s\n','Modulus','Real','Imaginary'))
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 0c115e4fecc528ca070168be0dd4947dc18aabf3..e254f4d85ae828e385ca25c2d375ee85ba7c3123 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -44,8 +44,8 @@ options_.gstep = ones(2,1);
 options_.gstep(1) = 1e-2;
 options_.gstep(2) = 1.0;
 options_.scalv = 1;
-options_.debug = 0;
-options_.initval_file = 0;
+options_.debug = false;
+options_.initval_file = false;
 options_.Schur_vec_tol = 1e-11; % used to find nonstationary variables in Schur decomposition of the
                                 % transition matrix
 options_.qz_criterium = [];
@@ -54,17 +54,17 @@ options_.lyapunov_complex_threshold = 1e-15;
 options_.solve_tolf = eps^(1/3);
 options_.solve_tolx = eps^(2/3);
 options_.dr_display_tol=1e-6;
-options_.minimal_workspace = 0;
+options_.minimal_workspace = false;
 options_.dp.maxit = 3000;
 options_.steady.maxit = 50;
 options_.simul.maxit = 50;
-options_.simul.robust_lin_solve = 0;
+options_.simul.robust_lin_solve = false;
 
-options_.mode_check.status = 0;
+options_.mode_check.status = false;
 options_.mode_check.neighbourhood_size = .5;
-options_.mode_check.symmetric_plots = 1;
+options_.mode_check.symmetric_plots = true;
 options_.mode_check.number_of_points = 20;
-options_.mode_check.nolik = 0;
+options_.mode_check.nolik = false;
 
 options_.huge_number = 1e7;
 
@@ -74,7 +74,7 @@ options_.threads.kronecker.sparse_hessian_times_B_kronecker_C = 1;
 options_.threads.local_state_space_iteration_2 = 1;
 
 % steady state
-options_.jacobian_flag = 1;
+options_.jacobian_flag = true;
 
 % steady state file
 if exist(['+' M_.fname '/steadystate.m'],'file')
@@ -85,7 +85,7 @@ else
     options_.steadystate_flag = 0;
 end
 options_.steadystate_partial = [];
-options_.steadystate.nocheck = 0;
+options_.steadystate.nocheck = false;
 
 % subset of the estimated deep parameters
 options_.ParamSubSet = 'None';
@@ -101,7 +101,7 @@ options_.bvar_prior_decay = 0.5;
 options_.bvar_prior_lambda = 5;
 options_.bvar_prior_mu = 2;
 options_.bvar_prior_omega = 1;
-options_.bvar_prior_flat = 0;
+options_.bvar_prior_flat = false;
 options_.bvar_prior_train = 0;
 options_.bvar.conf_sig = 0.6;
 
@@ -118,7 +118,7 @@ gmhmaxlik.target = 1/3; % Target for the acceptance rate.
 options_.gmhmaxlik = gmhmaxlik;
 
 % Request user input.
-options_.nointeractive = 0;
+options_.nointeractive = false;
 
 % Graphics
 options_.graphics.nrows = 3;
@@ -126,50 +126,50 @@ options_.graphics.ncols = 3;
 options_.graphics.line_types = {'b-'};
 options_.graphics.line_width = 1;
 options_.graph_format = 'eps';
-options_.nodisplay = 0;
-options_.nograph = 0;
-options_.no_graph.posterior = 0;
-options_.no_graph.shock_decomposition = 0;
+options_.nodisplay = false;
+options_.nograph = false;
+options_.no_graph.posterior = false;
+options_.no_graph.shock_decomposition = false;
 options_.XTick = [];
 options_.XTickLabel = [];
-options_.console_mode = 0;
+options_.console_mode = false;
 if isoctave
     if sum(get(0,'screensize'))==4
-        options_.console_mode = 1;
-        options_.nodisplay = 1;
+        options_.console_mode = true;
+        options_.nodisplay = true;
     end
 else
     if isunix && (~usejava('jvm') || ~feature('ShowFigureWindows'))
-        options_.console_mode = 1;
-        options_.nodisplay = 1;
+        options_.console_mode = true;
+        options_.nodisplay = true;
     end
 end
 
 % IRFs & other stoch_simul output
 options_.irf = 40;
 options_.impulse_responses.plot_threshold=1e-10;
-options_.relative_irf = 0;
+options_.relative_irf = false;
 options_.ar = 5;
 options_.hp_filter = 0;
 options_.one_sided_hp_filter = 0;
 options_.hp_ngrid = 512;
-options_.nodecomposition = 0;
-options_.nomoments = 0;
-options_.nocorr = 0;
+options_.nodecomposition = false;
+options_.nomoments = false;
+options_.nocorr = false;
 options_.periods = 0;
-options_.noprint = 0;
-options_.SpectralDensity.trigger = 0;
+options_.noprint = false;
+options_.SpectralDensity.trigger = false;
 options_.SpectralDensity.plot  = 1;
 options_.SpectralDensity.cutoff  = 150;
 options_.SpectralDensity.sdl = 0.01;
 options_.nofunctions = false;
 
-options_.bandpass.indicator = 0;
+options_.bandpass.indicator = false;
 options_.bandpass.passband = [6; 32];
 options_.bandpass.K=12;
 
-options_.irf_opt.diagonal_only = 0;
-options_.irf_opt.stderr_multiples = 0;
+options_.irf_opt.diagonal_only = false;
+options_.irf_opt.stderr_multiples = false;
 options_.irf_opt.irf_shock_graphtitles = {};
 options_.irf_opt.irf_shocks = [];
 
@@ -237,7 +237,7 @@ options_.bnlms = bnlms;
 % Particle filter
 %
 % Default is that we do not use the non linear kalman filter
-particle.status = 0;
+particle.status = false;
 % How do we initialize the states?
 particle.initialization = 1;
 particle.initial_state_prior_std = .1;
@@ -256,23 +256,23 @@ particle.unscented.alpha = 1;
 particle.unscented.beta = 2;
 particle.unscented.kappa = 1;
 % Configuration of resampling in case of particles
-particle.resampling.status.systematic = 1;
-particle.resampling.status.none = 0;
-particle.resampling.status.generic = 0;
+particle.resampling.status.systematic = true;
+particle.resampling.status.none = false;
+particle.resampling.status.generic = false;
 particle.resampling.threshold = .5;
-particle.resampling.method.kitagawa = 1;
-particle.resampling.method.smooth = 0;
-particle.resampling.method.stratified = 0;
+particle.resampling.method.kitagawa = true;
+particle.resampling.method.smooth = false;
+particle.resampling.method.stratified = false;
 % Set default algorithm
 particle.filter_algorithm = 'sis';
 % Approximation of the proposal distribution
-particle.proposal_approximation.cubature = 0;
-particle.proposal_approximation.unscented = 1;
-particle.proposal_approximation.montecarlo = 0;
+particle.proposal_approximation.cubature = false;
+particle.proposal_approximation.unscented = true;
+particle.proposal_approximation.montecarlo = false;
 % Approximation of the particle distribution
-particle.distribution_approximation.cubature = 0;
-particle.distribution_approximation.unscented = 1;
-particle.distribution_approximation.montecarlo = 0;
+particle.distribution_approximation.cubature = false;
+particle.distribution_approximation.unscented = true;
+particle.distribution_approximation.montecarlo = false;
 % Number of partitions for the smoothed resampling method
 particle.resampling.number_of_partitions = 200;
 % Configuration of the mixture filters
@@ -285,8 +285,8 @@ particle.mixture_measurement_shocks = 1 ;
 particle.liu_west_delta = 0.99 ;
 particle.liu_west_chol_sigma_bar = .01 ;
 % Options for setting the weights in conditional particle filters.
-particle.cpf_weights_method.amisanotristani = 1;
-particle.cpf_weights_method.murrayjonesparslow = 0;
+particle.cpf_weights_method.amisanotristani = true;
+particle.cpf_weights_method.murrayjonesparslow = false;
 % Copy ep structure in options_ global structure
 options_.particle = particle;
 options_.rwgmh.init_scale = 1e-4 ;
@@ -294,7 +294,7 @@ options_.rwgmh.scale_chain = 1 ;
 options_.rwgmh.scale_shock = 1e-5 ;
 
 % TeX output
-options_.TeX = 0;
+options_.TeX = false;
 
 % Exel
 options_.xls_sheet = 1; % Octave does not support the empty string, rather use first sheet
@@ -311,31 +311,30 @@ options_.forecasts.conf_sig = 0.9;
 options_.conditional_forecast.conf_sig = 0.9;
 
 % Model
-options_.linear = 0;
+options_.linear = false;
 
 % Deterministic simulation
 options_.stack_solve_algo = 0;
 options_.markowitz = 0.5;
 options_.minimal_solving_periods = 1;
-options_.endogenous_terminal_period = 0;
-options_.no_homotopy = 0;
+options_.endogenous_terminal_period = false;
+options_.no_homotopy = false;
 
 % Solution
 options_.order = 2;
-options_.pruning = 0;
+options_.pruning = false;
 options_.solve_algo = 4;
-options_.linear = 0;
 options_.replic = 50;
 options_.simul_replic = 1;
 options_.drop = 100;
-options_.aim_solver = 0; % i.e. by default do not use G.Anderson's AIM solver, use mjdgges instead
-options_.k_order_solver=0; % by default do not use k_order_perturbation but mjdgges
-options_.partial_information = 0;
-options_.ACES_solver = 0;
+options_.aim_solver = false; % i.e. by default do not use G.Anderson's AIM solver, use mjdgges instead
+options_.k_order_solver = false; % by default do not use k_order_perturbation but mjdgges
+options_.partial_information = false;
+options_.ACES_solver = false;
 options_.conditional_variance_decomposition = [];
 
 % Ramsey policy
-options_.ramsey_policy = 0;
+options_.ramsey_policy = false;
 options_.instruments = {};
 options_.timeless = 0;
 options_.ramsey.maxit = 500;
@@ -352,30 +351,30 @@ options_.dataset.xls_range = [];
 options_.Harvey_scale_factor = 10;
 options_.MaxNumberOfBytes = 1e8;
 options_.MaximumNumberOfMegaBytes = 111;
-options_.analytic_derivation = 0;
+options_.analytic_derivation = 0; % Not a boolean, can also take values -1 or 2
 options_.analytic_derivation_mode = 0;
-options_.bayesian_irf = 0;
+options_.bayesian_irf = false;
 options_.bayesian_th_moments = 0;
-options_.diffuse_filter = 0;
+options_.diffuse_filter = false;
 options_.filter_step_ahead = [];
-options_.filtered_vars = 0;
-options_.smoothed_state_uncertainty = 0;
+options_.filtered_vars = false;
+options_.smoothed_state_uncertainty = false;
 options_.first_obs = NaN;
 options_.nobs = NaN;
 options_.kalman_algo = 0;
-options_.fast_kalman_filter = 0;
+options_.fast_kalman_filter = false;
 options_.kalman_tol = 1e-10;
-options_.kalman.keep_kalman_algo_if_singularity_is_detected = 0;
+options_.kalman.keep_kalman_algo_if_singularity_is_detected = false;
 options_.diffuse_kalman_tol = 1e-6;
 options_.use_univariate_filters_if_singularity_is_detected = 1;
 options_.riccati_tol = 1e-6;
 options_.lik_algo = 1;
 options_.lik_init = 1;
-options_.load_mh_file = 0;
-options_.load_results_after_load_mh = 0;
-options_.logdata = 0;
-options_.loglinear = 0;
-options_.linear_approximation = 0;
+options_.load_mh_file = false;
+options_.load_results_after_load_mh = false;
+options_.logdata = false;
+options_.loglinear = false;
+options_.linear_approximation = false;
 options_.logged_steady_state = 0;
 options_.mh_conf_sig = 0.90;
 options_.prior_interval = 0.90;
@@ -393,7 +392,7 @@ options_.mh_tune_jscale.c3 = 4;
 options_.mh_init_scale = 2*options_.mh_jscale;
 options_.mh_mode = 1;
 options_.mh_nblck = 2;
-options_.mh_recover = 0;
+options_.mh_recover = false;
 options_.mh_replic = 20000;
 options_.recursive_estimation_restart = 0;
 options_.MCMC_jumping_covariance='hessian';
@@ -401,7 +400,7 @@ options_.use_calibration_initialization = 0;
 options_.endo_vars_for_moment_computations_in_estimation=[];
 
 % Run optimizer silently
-options_.silent_optimizer=0;
+options_.silent_optimizer = false;
 
 % Prior restrictions
 options_.prior_restrictions.status = 0;
@@ -409,15 +408,15 @@ options_.prior_restrictions.routine = [];
 
 options_.mode_compute = 4;
 options_.mode_file = '';
-options_.moments_varendo = 0;
+options_.moments_varendo = false;
 options_.nk = 1;
-options_.noconstant = 0;
-options_.nodiagnostic = 0;
+options_.noconstant = false;
+options_.nodiagnostic = false;
 options_.mh_posterior_mode_estimation = 0;
 options_.prefilter = 0;
 options_.presample = 0;
 options_.prior_trunc = 1e-10;
-options_.smoother = 0;
+options_.smoother = false;
 options_.posterior_max_subsample_draws = 1200;
 options_.sub_draws = [];
 options_.ME_plot_tol=1e-6;
@@ -469,13 +468,13 @@ for i=1:length(years)
     options_.conditional_variance_decomposition_dates(i) = ...
         (years(i)-1)*4+quarter;
 end
-options_.filter_covariance = 0;
-options_.filter_decomposition = 0;
-options_.selected_variables_only = 0;
-options_.contemporaneous_correlation = 0;
+options_.filter_covariance = false;
+options_.filter_decomposition = false;
+options_.selected_variables_only = false;
+options_.contemporaneous_correlation = false;
 options_.initialize_estimated_parameters_with_the_prior_mode = 0;
-options_.estimation_dll = 0;
-options_.estimation.moments_posterior_density.indicator = 1;
+options_.estimation_dll = false;
+options_.estimation.moments_posterior_density.indicator = true;
 options_.estimation.moments_posterior_density.gridpoints = 2^9;
 options_.estimation.moments_posterior_density.bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
 options_.estimation.moments_posterior_density.kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
@@ -485,7 +484,7 @@ options_.estimation.moments_posterior_density.kernel_function = 'gaussian'; % Ga
 % homotopy for steady state
 options_.homotopy_mode = 0;
 options_.homotopy_steps = 1;
-options_.homotopy_force_continue = 0;
+options_.homotopy_force_continue = false;
 
 % numerical hessian
 hessian.use_penalized_objective = false;
@@ -597,39 +596,39 @@ options_.prior_mc = 20000;
 options_.prior_analysis_endo_var_list = {};
 
 % did model undergo block decomposition + minimum feedback set computation ?
-options_.block = 0;
+options_.block = false;
 
 % model evaluated using a compiled MEX
-options_.use_dll = 0;
+options_.use_dll = false;
 
 % model evaluated using bytecode.dll
-options_.bytecode = 0;
+options_.bytecode = false;
 
-% if equal to 1 use a fixed point method to solve Sylvester equation (for large scale models)
-options_.sylvester_fp = 0;
+% if true, use a fixed point method to solve Sylvester equation (for large scale models)
+options_.sylvester_fp = false;
 
 % convergence criteria to solve iteratively a sylvester equations
 options_.sylvester_fixed_point_tol = 1e-12;
 
-% if 1 use a fixed point method to solve Lyapunov equation (for large scale models)
-options_.lyapunov_fp = 0;
-% if 1 use a doubling algorithm to solve Lyapunov equation (for large scale models)
-options_.lyapunov_db = 0;
-% if 1 use a square root solver to solve Lyapunov equation (for large scale models)
-options_.lyapunov_srs = 0;
+% if true, use a fixed point method to solve Lyapunov equation (for large scale models)
+options_.lyapunov_fp = false;
+% if true, use a doubling algorithm to solve Lyapunov equation (for large scale models)
+options_.lyapunov_db = false;
+% if true, use a square root solver to solve Lyapunov equation (for large scale models)
+options_.lyapunov_srs = false;
 
 % convergence criterion for iteratives methods to solve lyapunov equations
 options_.lyapunov_fixed_point_tol = 1e-10;
 options_.lyapunov_doubling_tol = 1e-16;
 
-% if equal to 1 use a cycle reduction method to compute the decision rule (for large scale models)
-options_.dr_cycle_reduction = 0;
+% if true, use a cycle reduction method to compute the decision rule (for large scale models)
+options_.dr_cycle_reduction = false;
 
 % convergence criterion for iteratives methods to solve the decision rule
 options_.dr_cycle_reduction_tol = 1e-7;
 
-% if equal to 1 use a logarithmic reduction method to compute the decision rule (for large scale models)
-options_.dr_logarithmic_reduction = 0;
+% if true, use a logarithmic reduction method to compute the decision rule (for large scale models)
+options_.dr_logarithmic_reduction = false;
 
 % convergence criterion for iteratives methods to solve the decision rule
 options_.dr_logarithmic_reduction_tol = 1e-12;
@@ -663,8 +662,8 @@ options_.nonlinear_filter = [];
 % SBVAR
 options_.ms.vlistlog = [];
 options_.ms.restriction_fname = 0;
-options_.ms.cross_restrictions = 0;
-options_.ms.contemp_reduced_form = 0;
+options_.ms.cross_restrictions = false;
+options_.ms.contemp_reduced_form = false;
 options_.ms.real_pseudo_forecast = 0;
 options_.ms.dummy_obs = 0;
 options_.ms.ncsk = 0;
@@ -695,7 +694,7 @@ options_.graph_save_formats.fig = 0;
 options_.risky_steadystate = 0;
 
 % endogenous prior
-options_.endogenous_prior = 0;
+options_.endogenous_prior = false;
 options_.endogenous_prior_restrictions.irf={};
 options_.endogenous_prior_restrictions.moment={};
 
@@ -703,17 +702,17 @@ options_.endogenous_prior_restrictions.moment={};
 options_.osr.opt_algo=4;
 
 % use GPU
-options_.gpu = 0;
+options_.gpu = false;
 
 %Geweke convergence diagnostics
 options_.convergence.geweke.taper_steps=[4 8 15];
 options_.convergence.geweke.geweke_interval=[0.2 0.5];
 %Raftery/Lewis convergence diagnostics;
-options_.convergence.rafterylewis.indicator=0;
+options_.convergence.rafterylewis.indicator=false;
 options_.convergence.rafterylewis.qrs=[0.025 0.005 0.95];
 
 % Options for lmmcp solver
-options_.lmmcp.status = 0;
+options_.lmmcp.status = false;
 
 % Options for lcppath solver
 options_.lcppath.A = [];
diff --git a/matlab/discretionary_policy.m b/matlab/discretionary_policy.m
index 6513cb06bb35f48c6a0d0f689347333adc93c076..d14f4044a47cc6cd36eb372d4cf0ee61598d9242 100644
--- a/matlab/discretionary_policy.m
+++ b/matlab/discretionary_policy.m
@@ -24,7 +24,7 @@ oldoptions = options_;
 options_.order = 1;
 info = stoch_simul(var_list);
 
-if options_.noprint == 0
+if ~options_.noprint
     disp_steady_state(M_,oo_)
     for i=M_.orig_endo_nbr:M_.endo_nbr
         if strmatch('mult_', M_.endo_names{i})
@@ -36,4 +36,4 @@ end
 
 oo_.planner_objective_value = evaluate_planner_objective(M_,options_,oo_);
 
-options_ = oldoptions;
\ No newline at end of file
+options_ = oldoptions;
diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m
index 938f0cdb2405629ba496c7ade16e5bc8bb1188b0..bd9744d362d0ad9bdee1477f068f6e916d7f344e 100644
--- a/matlab/disp_moments.m
+++ b/matlab/disp_moments.m
@@ -86,7 +86,7 @@ oo_.kurtosis = (mean(y.^4)./(s2.*s2)-3)';
 labels = M_.endo_names(ivar);
 labels_TeX = M_.endo_names_tex(ivar);
 
-if options_.nomoments == 0
+if ~options_.nomoments
     z = [ m' s' s2' (mean(y.^3)./s2.^1.5)' (mean(y.^4)./(s2.*s2)-3)' ];
     title='MOMENTS OF SIMULATED VARIABLES';
     title=add_filter_subtitle(title, options_);
@@ -97,12 +97,12 @@ if options_.nomoments == 0
     end
 end
 
-if options_.nocorr == 0
+if ~options_.nocorr
     corr = (y'*y/size(y,1))./(s'*s);
     if options_.contemporaneous_correlation
         oo_.contemporaneous_correlation = corr;
     end
-    if options_.noprint == 0
+    if ~options_.noprint
         title = 'CORRELATION OF SIMULATED VARIABLES';
         title=add_filter_subtitle(title,options_);
         headers = vertcat('VARIABLE', M_.endo_names(ivar));
@@ -115,7 +115,7 @@ if options_.nocorr == 0
     end
 end
 
-if options_.noprint == 0 && length(options_.conditional_variance_decomposition)
+if ~options_.noprint && length(options_.conditional_variance_decomposition)
     fprintf('\nSTOCH_SIMUL: conditional_variance_decomposition requires theoretical moments, i.e. periods=0.\n')
 end
 
@@ -126,7 +126,7 @@ if ar > 0
         oo_.autocorr{i} = y(ar+1:end,:)'*y(ar+1-i:end-i,:)./((size(y,1)-ar)*std(y(ar+1:end,:))'*std(y(ar+1-i:end-i,:)));
         autocorr = [ autocorr diag(oo_.autocorr{i}) ];
     end
-    if options_.noprint == 0
+    if ~options_.noprint
         title = 'AUTOCORRELATION OF SIMULATED VARIABLES';
         title=add_filter_subtitle(title,options_);
         headers = vertcat('VARIABLE', cellstr(int2str([1:ar]')));
@@ -236,4 +236,4 @@ else
     error('disp_moments:: You cannot use more than one filter at the same time')
 end
 
-end
\ No newline at end of file
+end
diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m
index ec9f6028be706ae5f0f217e9d9ec2838bae46055..fd99502d9763c1b1c9614e09809e333611ce65a9 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/disp_th_moments.m
@@ -141,7 +141,7 @@ if size(stationary_vars, 1) > 0
         StateSpaceModel.observable_pos = options_.varobs_id;
         [oo_.conditional_variance_decomposition, oo_.conditional_variance_decomposition_ME] = ...
             conditional_variance_decomposition(StateSpaceModel, conditional_variance_steps, ivar);
-        if options_.noprint == 0
+        if ~options_.noprint
             display_conditional_variance_decomposition(oo_.conditional_variance_decomposition, conditional_variance_steps, ivar, M_, options_);
             if ME_present
                 display_conditional_variance_decomposition(oo_.conditional_variance_decomposition_ME, conditional_variance_steps, ...
@@ -158,7 +158,7 @@ if length(i1) == 0
     return
 end
 
-if options_.nocorr == 0 && size(stationary_vars, 1)>0
+if ~options_.nocorr && size(stationary_vars, 1)>0
     corr = NaN(size(oo_.gamma_y{1}));
     corr(i1,i1) = oo_.gamma_y{1}(i1,i1)./(sd(i1)*sd(i1)');
     if options_.contemporaneous_correlation 
diff --git a/matlab/dr_block.m b/matlab/dr_block.m
index c0209d527fd892a567cdceeec8b765729cd4ab0f..0328c1048c24ec3432562717e137cc332a2eed34 100644
--- a/matlab/dr_block.m
+++ b/matlab/dr_block.m
@@ -425,7 +425,7 @@ for i = 1:Size
 
         row_indx = n_static+1:n;
 
-        if task ~= 1 && options_.dr_cycle_reduction == 1
+        if task ~= 1 && options_.dr_cycle_reduction
             A1 = [aa(row_indx,index_m ) zeros(n_dynamic,n_fwrd)];
             B1 = [aa(row_indx,index_0m) aa(row_indx,index_0p) ];
             C1 = [zeros(n_dynamic,n_pred) aa(row_indx,index_p)];
@@ -436,7 +436,7 @@ for i = 1:Size
             gx = ghx(1+n_pred:end,:);
         end
 
-        if (task ~= 1 && ((options_.dr_cycle_reduction == 1 && info ==1) || options_.dr_cycle_reduction == 0)) || task == 1
+        if (task ~= 1 && ((options_.dr_cycle_reduction && info ==1) || ~options_.dr_cycle_reduction)) || task == 1
             D = [[aa(row_indx,index_0m) zeros(n_dynamic,n_both) aa(row_indx,index_p)] ; [zeros(n_both, n_pred) eye(n_both) zeros(n_both, n_both + n_fwrd)]];
             E = [-aa(row_indx,[index_m index_0p])  ; [zeros(n_both, n_both + n_pred) eye(n_both, n_both + n_fwrd) ] ];
 
@@ -588,7 +588,7 @@ for i = 1:Size
                 if block_type == 5
                     vghx_other = - inv(kron(eye(size(D_,2)), A_) + kron(C_', B_)) * vec(D_);
                     ghx_other = reshape(vghx_other, size(D_,1), size(D_,2));
-                elseif options_.sylvester_fp == 1
+                elseif options_.sylvester_fp
                     ghx_other = gensylv_fp(A_, B_, C_, D_, i, options_.sylvester_fixed_point_tol);
                 else
                     [err, ghx_other] = gensylv(1, A_, B_, C_, -D_);
diff --git a/matlab/dyn_first_order_solver.m b/matlab/dyn_first_order_solver.m
index 861356a6321c820e94c928af260a9c97e5cb304e..63cb4f9da5a00fd13edcb9935b76f74a3d819b3b 100644
--- a/matlab/dyn_first_order_solver.m
+++ b/matlab/dyn_first_order_solver.m
@@ -167,7 +167,7 @@ if task ~= 1 && (DynareOptions.dr_cycle_reduction || DynareOptions.dr_logarithmi
     A1 = [aa(row_indx,index_m ) zeros(ndynamic,nfwrd)];
     B1 = [aa(row_indx,index_0m) aa(row_indx,index_0p) ];
     C1 = [zeros(ndynamic,npred) aa(row_indx,index_p)];
-    if DynareOptions.dr_cycle_reduction == 1
+    if DynareOptions.dr_cycle_reduction
         [ghx, info] = cycle_reduction(A1, B1, C1, DynareOptions.dr_cycle_reduction_tol);
     else
         [ghx, info] = logarithmic_reduction(C1, B1, A1, DynareOptions.dr_logarithmic_reduction_tol, DynareOptions.dr_logarithmic_reduction_maxiter);
diff --git a/matlab/dyn_forecast.m b/matlab/dyn_forecast.m
index 162ace08ca02043562920cf99824893d0698f966..f152c879fa1ccadf708a4ff907cf656a89603925 100644
--- a/matlab/dyn_forecast.m
+++ b/matlab/dyn_forecast.m
@@ -169,7 +169,7 @@ if ~isscalar(trend) %add trend back to forecast
     yf(i_var_obs,:) = yf(i_var_obs,:) + trend;
 end
 
-if options.loglinear == 1
+if options.loglinear
     if options.prefilter == 1 %subtract steady state and add mean for observables
         yf(i_var_obs,:)=yf(i_var_obs,:)-repmat(log(oo.dr.ys(i_var_obs)),1,horizon+M.maximum_lag)+ repmat(mean_varobs,1,horizon+M.maximum_lag);
     end
@@ -194,7 +194,7 @@ for i=1:M.exo_det_nbr
     forecast.Exogenous.(M.exo_det_names{i}) = oo.exo_det_simul(maximum_lag+(1:horizon),i);
 end
 
-if options.nograph == 0
+if ~options.nograph
     oo.forecast = forecast;
     forecast_graphs(var_list, M, oo, options)
 end
diff --git a/matlab/dyn_ramsey_static.m b/matlab/dyn_ramsey_static.m
index ccbad5e7bc57108297f292e5d7bd14d8ca65c7ad..e94deb5044772f6df93fd26f244be6269ed430e0 100644
--- a/matlab/dyn_ramsey_static.m
+++ b/matlab/dyn_ramsey_static.m
@@ -65,7 +65,7 @@ elseif options_.steadystate_flag
         %solve for instrument, using multivariate solver, starting at
         %initial value for instrument
         opt = options_;
-        opt.jacobian_flag = 0;
+        opt.jacobian_flag = false;
         [inst_val,info1] = dynare_solve(nl_func,ys_init(k_inst), ...
                                         opt);
         if info1~=0
@@ -80,7 +80,7 @@ else
     n_var = M.orig_endo_nbr;
     xx = oo.steady_state(1:n_var);
     opt = options_;
-    opt.jacobian_flag = 0;
+    opt.jacobian_flag = false;
     [xx,info1] = dynare_solve(nl_func,xx,opt);
     if info1~=0
         check=81;
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 776adbaa05953f0347b303870037faf004daed67..933321d9bb2de7f1cbc90d19bd74b39f149825c7 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -81,7 +81,7 @@ if options_.order > 1
     elseif options_.particle.status && options_.order>2
         error(['Non linear filter are not implemented with order ' int2str(options_.order) ' approximation of the model!'])
     elseif ~options_.particle.status && options_.order==2
-        error('For estimating the model with a second order approximation using a non linear filter, one should have options_.particle.status=1;')
+        error('For estimating the model with a second order approximation using a non linear filter, one should have options_.particle.status=true;')
     else
         error(['Cannot estimate a model with an order ' int2str(options_.order) ' approximation!'])
     end
@@ -129,9 +129,9 @@ if options_.dsge_var
 end
 
 % Set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file).
-M_.sigma_e_is_diagonal = 1;
+M_.sigma_e_is_diagonal = true;
 if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(estim_params_,'calibrated_covariances')
-    M_.sigma_e_is_diagonal = 0;
+    M_.sigma_e_is_diagonal = false;
 end
 
 data = dataset_.data;
@@ -181,7 +181,7 @@ catch % if check fails, provide info on using calibration if present
 end
 
 if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
-    if options_.smoother == 1
+    if options_.smoother
         [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
         [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
     end
@@ -320,7 +320,7 @@ if ~options_.mh_posterior_mode_estimation && options_.cova_compute
     end
 end
 
-if options_.mode_check.status == 1 && ~options_.mh_posterior_mode_estimation
+if options_.mode_check.status && ~options_.mh_posterior_mode_estimation
     ana_deriv_old = options_.analytic_derivation;
     options_.analytic_derivation = 0;
     mode_check(objective_function,xparam1,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
@@ -549,7 +549,7 @@ if options_.particle.status
 end
 
 if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file)) ...
-    || ~options_.smoother ) && options_.partial_information == 0  % to be fixed
+    || ~options_.smoother ) && ~options_.partial_information  % to be fixed
     %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothes variables
     [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_);
     [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
@@ -799,4 +799,4 @@ if reset_options_related_to_estimation
 end
 if first_obs_nan_indicator
     options_.first_obs=NaN;
-end
\ No newline at end of file
+end
diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index 7676c644575736644517fca349c0ba9fbe794e5b..5011a392aa484d837e27659e04b2e70479e4df67 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -448,7 +448,7 @@ else
 end
 
 % Define union of observed and state variables
-if options_.block == 1
+if options_.block
     k1 = k1';
     [k2, i_posA, i_posB] = union(k1', M_.state_var', 'rows');
     % Set restrict_state to postion of observed + state variables in expanded state vector.
diff --git a/matlab/dynare_resolve.m b/matlab/dynare_resolve.m
index 32a1a885d4497160e5b67fe0a51a2c3770b0209a..5fcbc7151c921ade1a26353d3ab272bfa0559dbc 100644
--- a/matlab/dynare_resolve.m
+++ b/matlab/dynare_resolve.m
@@ -86,7 +86,7 @@ switch nargin
     nstatic = Model.nstatic;
     nspred = Model.nspred;
     iv = (1:endo_nbr)';
-    if DynareOptions.block == 0
+    if ~DynareOptions.block
         ic = [ nstatic+(1:nspred) endo_nbr+(1:size(DynareResults.dr.ghx,2)-nspred) ]';
     else
         ic = DynareResults.dr.restrict_columns;
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 4dc408eb6ea8eed0472c3b6c4906f2a4057bb072..ebb883c9adbc253e3b5a4bf297d765a528015af0 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -34,8 +34,8 @@ function [x,info,fvec,fjac] = dynare_solve(func,x,options,varargin)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-% jacobian_flag=1:  jacobian given by the 'func' function
-% jacobian_flag=0:  jacobian obtained numerically
+% jacobian_flag=true:  jacobian given by the 'func' function
+% jacobian_flag=false:  jacobian obtained numerically
 jacobian_flag = options.jacobian_flag;
 
 % Set tolerance parameter depending the the caller function.
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index d8120d6610c23c6cb9b4d9f3ef08831b70e5b957..12507a2b4686e964a6707b1668e3316553015d68 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -214,8 +214,8 @@ elseif steadystate_flag
     if info(1)
         return
     end
-elseif (options.bytecode == 0 && options.block == 0)
-    if options.linear == 0
+elseif ~options.bytecode && ~options.block
+    if ~options.linear
         % non linear model
         static_model = str2func([M.fname '.static']);
         [ys,check] = dynare_solve(@static_problem,...
diff --git a/matlab/initvalf.m b/matlab/initvalf.m
index 33d0e3a022e0935b489b33e9d7482813594c4ada..f301e045fd5a8f927a7adfa62519abecfc034d15 100644
--- a/matlab/initvalf.m
+++ b/matlab/initvalf.m
@@ -69,7 +69,7 @@ switch (extension)
     error(['Unsupported extension for datafile: ' extension])
 end
 
-options_.initval_file = 1;
+options_.initval_file = true;
 oo_.endo_simul = [];
 oo_.exo_simul = [];
 
@@ -107,4 +107,4 @@ for i_=1:length(M_.exo_names)
         x_ = data_(:,k_);
         oo_.exo_simul = [oo_.exo_simul x_];
     end
-end
\ No newline at end of file
+end
diff --git a/matlab/lyapunov_solver.m b/matlab/lyapunov_solver.m
index 2bafb065ab159cefbfaa8c3a22666f0800ab69c9..3fa9dd42d95471b6b3442c4a0af1404b909df02d 100644
--- a/matlab/lyapunov_solver.m
+++ b/matlab/lyapunov_solver.m
@@ -15,11 +15,11 @@ function P=lyapunov_solver(T,R,Q,DynareOptions) % --*-- Unitary tests --*--
 % Algorithms
 %   Default, if none of the other algorithms is selected:
 %       Reordered Schur decomposition (Bartels-Stewart algorithm)
-%   DynareOptions.lyapunov_fp == 1
+%   DynareOptions.lyapunov_fp == true
 %       iteration-based fixed point algorithm
-%   DynareOptions.lyapunov_db == 1
+%   DynareOptions.lyapunov_db == true
 %       doubling algorithm
-%   DynareOptions.lyapunov_srs == 1
+%   DynareOptions.lyapunov_srs == true
 %       Square-root solver for discrete-time Lyapunov equations (requires Matlab System Control toolbox
 %       or Octave control package)
 
@@ -40,14 +40,14 @@ function P=lyapunov_solver(T,R,Q,DynareOptions) % --*-- Unitary tests --*--
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-if DynareOptions.lyapunov_fp == 1
+if DynareOptions.lyapunov_fp
     P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, DynareOptions.debug);
-elseif DynareOptions.lyapunov_db == 1
+elseif DynareOptions.lyapunov_db
     [P, errorflag] = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol);
     if errorflag %use Schur-based method
         P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug);
     end
-elseif DynareOptions.lyapunov_srs == 1
+elseif DynareOptions.lyapunov_srs
     % works only with Matlab System Control toolbox or Octave control package,
     if isoctave
         if ~user_has_octave_forge_package('control')
@@ -72,7 +72,7 @@ end
 %$ options_.qz_criterium=1-options_.qz_zero_threshold;
 %$ options_.lyapunov_fixed_point_tol = 1e-10;
 %$ options_.lyapunov_doubling_tol = 1e-16;
-%$ options_.debug=0;
+%$ options_.debug=false;
 %$
 %$ n_small=8;
 %$ m_small=10;
@@ -91,7 +91,7 @@ end
 %$ R_large=randn(n_large,m_large);
 %$
 %$ % DynareOptions.lyapunov_fp == 1
-%$ options_.lyapunov_fp = 1;
+%$ options_.lyapunov_fp = true;
 %$ try
 %$    Pstar1_small = lyapunov_solver(T_small,R_small,Q_small,options_);
 %$    Pstar1_large = lyapunov_solver(T_large,R_large,Q_large,options_);
@@ -101,8 +101,8 @@ end
 %$ end
 %$
 %$ % Dynareoptions.lyapunov_db == 1
-%$ options_.lyapunov_fp = 0;
-%$ options_.lyapunov_db = 1;
+%$ options_.lyapunov_fp = false;
+%$ options_.lyapunov_db = true;
 %$ try
 %$    Pstar2_small = lyapunov_solver(T_small,R_small,Q_small,options_);
 %$    Pstar2_large = lyapunov_solver(T_large,R_large,Q_large,options_);
@@ -113,8 +113,8 @@ end
 %$
 %$ % Dynareoptions.lyapunov_srs == 1
 %$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
-%$     options_.lyapunov_db = 0;
-%$     options_.lyapunov_srs = 1;
+%$     options_.lyapunov_db = false;
+%$     options_.lyapunov_srs = true;
 %$     try
 %$        Pstar3_small = lyapunov_solver(T_small,R_small,Q_small,options_);
 %$        Pstar3_large = lyapunov_solver(T_large,R_large,Q_large,options_);
@@ -127,7 +127,7 @@ end
 %$ end
 %$
 %$ % Standard
-%$     options_.lyapunov_srs = 0;
+%$     options_.lyapunov_srs = false;
 %$ try
 %$    Pstar4_small = lyapunov_solver(T_small,R_small,Q_small,options_);
 %$    Pstar4_large = lyapunov_solver(T_large,R_large,Q_large,options_);
diff --git a/matlab/mode_check.m b/matlab/mode_check.m
index b9a24e3153d152d1cba2c559e15f09cdd7d0d704..ca0bb96ab778c4b27feb554847ddbf4c957cb17a 100644
--- a/matlab/mode_check.m
+++ b/matlab/mode_check.m
@@ -86,7 +86,7 @@ end
 
 ll = DynareOptions.mode_check.neighbourhood_size;
 if isinf(ll)
-    DynareOptions.mode_check.symmetric_plots = 0;
+    DynareOptions.mode_check.symmetric_plots = false;
 end
 
 mcheck = struct('cross',struct(),'emode',struct());
@@ -149,7 +149,7 @@ for plt = 1:nbplt
         z1 = l1:((x(kk)-l1)/(DynareOptions.mode_check.number_of_points/2)):x(kk);
         z2 = x(kk):((l2-x(kk))/(DynareOptions.mode_check.number_of_points/2)):l2;
         z  = union(z1,z2);
-        if DynareOptions.mode_check.nolik==0
+        if ~DynareOptions.mode_check.nolik
             y = zeros(length(z),2);
             dy = priordens(xx,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
         end
@@ -164,7 +164,7 @@ for plt = 1:nbplt
                     fprintf('mode_check:: could not solve model for parameter %s at value %4.3f, error code: %u\n',name,z(i),info(1))
                 end
             end
-            if DynareOptions.mode_check.nolik==0
+            if ~DynareOptions.mode_check.nolik
                 lnprior = priordens(xx,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
                 y(i,2)  = (y(i,1)+lnprior-dy);
             end
@@ -189,7 +189,7 @@ for plt = 1:nbplt
         hold off
         drawnow
     end
-    if DynareOptions.mode_check.nolik==0
+    if ~DynareOptions.mode_check.nolik
         if isoctave
             axes('outerposition',[0.3 0.93 0.42 0.07],'box','on'),
         else
diff --git a/matlab/ms-sbvar/initialize_ms_sbvar_options.m b/matlab/ms-sbvar/initialize_ms_sbvar_options.m
index b618aef3cd7b41d18b49e757dc506e56e0010b26..ec91ad1e3009a87f52591879d60e3a989540f98f 100644
--- a/matlab/ms-sbvar/initialize_ms_sbvar_options.m
+++ b/matlab/ms-sbvar/initialize_ms_sbvar_options.m
@@ -36,7 +36,7 @@ options_.ms.mh_file = '';
 options_.ms.free_param_file = '';
 options_.ms.output_file_tag = '';
 options_.ms.simulation_file_tag = '';
-options_.ms.create_init = 1;
+options_.ms.create_init = true;
 % prepare ms sbvar & estimation
 options_.ms.coefficients_prior_hyperparameters = [1.0 1.0 0.1 1.2 1.0 1.0];
 options_.ms.freq = 4;
@@ -45,7 +45,7 @@ options_.ms.final_subperiod = '';
 options_.ms.nlags = 1;
 options_.ms.cross_restrictions = 0;
 options_.ms.contemp_reduced_form = 0;
-options_.ms.bayesian_prior = 1;
+options_.ms.bayesian_prior = true;
 options_.ms.alpha = 1;
 options_.ms.beta = 1;
 options_.ms.gsig2_lmdm = 50^2;
@@ -73,28 +73,28 @@ options_.ms.mh_replic = 10000; % default differs from Dan's code
 options_.ms.thinning_factor = 1;
 options_.ms.drop = 0.1*options_.ms.mh_replic*options_.ms.thinning_factor;
 options_.ms.adaptive_mh_draws = 30000;
-options_.ms.save_draws = 0;
+options_.ms.save_draws = false;
 % mdd
 options_.ms.proposal_draws = 100000;
-options_.ms.use_mean_center = 0;
+options_.ms.use_mean_center = false;
 options_.ms.proposal_type = 3;
 options_.ms.proposal_lower_bound = 0.1;
 options_.ms.proposal_upper_bound = 0.9;
 % probabilities
 options_.ms.filtered_probabilities = 0;
-options_.ms.real_time_smoothed_probabilities = 0;
+options_.ms.real_time_smoothed_probabilities = false;
 % irf
 options_.ms.horizon = 12;
-options_.ms.filtered_probabilities = 0;
+options_.ms.filtered_probabilities = false;
 options_.ms.percentiles = [.16 .5 .84];
-options_.ms.parameter_uncertainty = 0;
+options_.ms.parameter_uncertainty = false;
 options_.ms.shock_draws = 10000;
 options_.ms.shocks_per_parameter = 10;
 options_.ms.median = 0;
 options_.ms.regime = 0;
-options_.ms.regimes = 0;
+options_.ms.regimes = false;
 % forecast
 options_.ms.forecast_data_obs = 0;
 % variance decomposition
-options_.ms.error_bands = 1;
+options_.ms.error_bands = true;
 end
diff --git a/matlab/mult_elimination.m b/matlab/mult_elimination.m
index 98c17aa7ca96dc3b28611a3b6cd6bac943f69f2b..5ff1418c3e336148f648d3685d5d46dcf53c3191 100644
--- a/matlab/mult_elimination.m
+++ b/matlab/mult_elimination.m
@@ -101,7 +101,7 @@ dr.M4 = M4;
 
 nvar = length(varlist);
 
-if nvar > 0 && options_.noprint == 0
+if nvar > 0 && ~options_.noprint
     res_table = zeros(2*(nm_nbr+M_.exo_nbr),nvar);
     headers = {'Variables'};
     for i=1:length(varlist)
diff --git a/matlab/partial_information/PCL_resol.m b/matlab/partial_information/PCL_resol.m
index a9e7fc01b969f8f3b314a9826bf164964e8501fc..9ddda8c13c33c635dcf910c6d5cb8df6620f01d9 100644
--- a/matlab/partial_information/PCL_resol.m
+++ b/matlab/partial_information/PCL_resol.m
@@ -45,7 +45,7 @@ function [dr,info]=PCL_resol(ys,check_flag)
 global M_ options_ oo_
 global it_
 
-jacobian_flag = 0;
+jacobian_flag = false;
 
 info = 0;
 
@@ -84,13 +84,13 @@ if options_.steadystate_flag
 
 else
     % testing if ys isn't a steady state or if we aren't computing Ramsey policy
-    if  options_.ramsey_policy == 0
+    if ~options_.ramsey_policy
         if options_.linear == 0
             % nonlinear models
             if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; ...
                                     oo_.exo_det_steady_state], M_.params))) > options_.dynatol.f
                 opt = options_;
-                opt.jacobian_flag = 0;
+                opt.jacobian_flag = false;
                 [dr.ys,check1] = dynare_solve(fh,dr.ys,opt,...
                                               [oo_.exo_steady_state; ...
                                     oo_.exo_det_steady_state], M_.params);
@@ -126,7 +126,7 @@ if ~isreal(dr.ys)
 end
 
 dr.fbias = zeros(M_.endo_nbr,1);
-if( (options_.partial_information ==1) || (options_.ACES_solver==1))%&& (check_flag == 0)
+if( options_.partial_information || options_.ACES_solver)%&& (check_flag == 0)
     [dr,info,M_,options_,oo_] = dr1_PI(dr,check_flag,M_,options_,oo_);
 else
     [dr,info,M_,options_,oo_] = dr1(dr,check_flag,M_,options_,oo_);
diff --git a/matlab/partial_information/PI_gensys.m b/matlab/partial_information/PI_gensys.m
index 8e905b8ff18dbd27ac30cac2fad5fcc8b6692e20..92d4831667528ab19592335f867d35907cad3556 100644
--- a/matlab/partial_information/PI_gensys.m
+++ b/matlab/partial_information/PI_gensys.m
@@ -89,7 +89,7 @@ try
     warning('on','MATLAB:singularMatrix');
     warning('on','MATLAB:nearlySingularMatrix');
     if (any(any(isinf(UAVinv))) || any(any(isnan(UAVinv))))
-        if(options_.ACES_solver==1)
+        if(options_.ACES_solver)
             disp('ERROR! saving PI_gensys_data_dump');
             save PI_gensys_data_dump
             error('PI_gensys: Inversion of poss. zero matrix UAVinv=inv(U02''*a1*V02)!');
@@ -146,7 +146,7 @@ G1pi=[Ze11 Ze12 Ze134 Ze134; P1 G11 G12 G13; Ze31 G21 G22 G23; P3 G31 G32 G33];
 
 impact=[eye(NX,NX); zeros(n+FL_RANK,NX)];
 
-if(options_.ACES_solver==1)
+if(options_.ACES_solver)
     if isfield(lq_instruments,'names')
         num_inst=size(lq_instruments.names,1);
         if num_inst>0
@@ -251,7 +251,7 @@ if zxz
     nmat=[]; %;gev=[]
     return
 end
-if (FL_RANK ~= nunstab && options_.ACES_solver~=1)
+if (FL_RANK ~= nunstab && ~options_.ACES_solver)
     disp(['Number of unstable variables ' num2str(nunstab)]);
     disp( ['does not match number of expectational equations ' num2str(FL_RANK)]);
     nmat=[];% gev=[];
diff --git a/matlab/partial_information/dr1_PI.m b/matlab/partial_information/dr1_PI.m
index b741854968b4c091a57ee05c591756a18a6b1819..2a65722684c3612dc51c94edd519ec2e157521b0 100644
--- a/matlab/partial_information/dr1_PI.m
+++ b/matlab/partial_information/dr1_PI.m
@@ -64,8 +64,8 @@ options_ = set_default_option(options_,'qz_criterium',1.000001);
 
 xlen = M_.maximum_endo_lead + M_.maximum_endo_lag + 1;
 
-if (options_.aim_solver == 1)
-    options_.aim_solver == 0;
+if options_.aim_solver
+    options_.aim_solver = false;
     warning('You can not use AIM with Part Info solver. AIM ignored');
 end
 if (options_.order > 1)
@@ -74,7 +74,7 @@ if (options_.order > 1)
 end
 
 % expanding system for Optimal Linear Regulator
-if options_.ramsey_policy && options_.ACES_solver == 0
+if options_.ramsey_policy && ~options_.ACES_solver
     if isfield(M_,'orig_model')
         orig_model = M_.orig_model;
         M_.endo_nbr = orig_model.endo_nbr;
@@ -88,7 +88,7 @@ if options_.ramsey_policy && options_.ACES_solver == 0
     old_solve_algo = options_.solve_algo;
     %  options_.solve_algo = 1;
     opt = options_;
-    opt.jacobian_flag = 0;
+    opt.jacobian_flag = false;
     oo_.steady_state = dynare_solve('ramsey_static',oo_.steady_state,opt,M_,options_,oo_,it_);
     options_.solve_algo = old_solve_algo;
     [~,~,multbar] = ramsey_static(oo_.steady_state,M_,options_,oo_,it_);
@@ -107,7 +107,7 @@ else
     end
 
 
-    if options_.ACES_solver == 1
+    if options_.ACES_solver
         sim_ruleids=[];
         tct_ruleids=[];
         if  size(M_.equations_tags,1)>0  % there are tagged equations, check if they are aceslq rules
@@ -133,7 +133,7 @@ else
     [~,jacobia_] = feval([M_.fname '.dynamic'],z,[oo_.exo_simul ...
                         oo_.exo_det_simul], M_.params, dr.ys, it_);
 
-    if options_.ACES_solver==1 && (length(sim_ruleids)>0 || length(tct_ruleids)>0 )
+    if options_.ACES_solver && (length(sim_ruleids)>0 || length(tct_ruleids)>0 )
         if length(sim_ruleids)>0
             sim_rule=jacobia_(sim_ruleids,:);
             % uses the subdirectory - BY
@@ -172,7 +172,7 @@ sdyn = M_.endo_nbr - nstatic;
 k0 = M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var);
 k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
 
-if (options_.aim_solver == 1)
+if options_.aim_solver
     error('Anderson and Moore AIM solver is not compatible with Partial Information models');
 end % end if useAIM and...
 
@@ -182,7 +182,7 @@ end % end if useAIM and...
     nendo=M_.endo_nbr; % = size(aa,1)
 
 
-    if(options_.ACES_solver==1)
+    if(options_.ACES_solver)
         %if ~isfield(lq_instruments,'names')
         if isfield(options_,'instruments')
             lq_instruments.names=options_.instruments;
@@ -262,7 +262,7 @@ end % end if useAIM and...
                 fnd = find(lead_lag(:,1));
                 AA2(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,1))); %backward
             end
-        elseif options_.ACES_solver==1 % more endo vars than equations in jacobia_
+        elseif options_.ACES_solver % more endo vars than equations in jacobia_
         if nendo-xlen==num_inst
             PSI=[PSI;zeros(num_inst, M_.exo_nbr)];
             % AA1 contemporary
@@ -324,7 +324,7 @@ end % end if useAIM and...
         %      (b) matrices TT1, TT2  that relate y(t) to these states:
         %      y(t)=[TT1 TT2][s(t)' x(t)']'.
 
-        if(options_.ACES_solver==1)
+        if(options_.ACES_solver)
             if isfield(lq_instruments,'xsopt_SS')
                 SSbar= diag([lq_instruments.xsopt_SS(m_var)]);% lq_instruments.xsopt_SS(lq_instruments.inst_var_indices)]);
                 insSSbar=repmat(lq_instruments.xsopt_SS(lq_instruments.inst_var_indices)',nendo-num_inst,1);
@@ -349,7 +349,7 @@ end % end if useAIM and...
         end
 
         % reuse some of the bypassed code and tests that may be needed
-        if (eu(1) ~= 1 || eu(2) ~= 1) && options_.ACES_solver==0
+        if (eu(1) ~= 1 || eu(2) ~= 1) && ~options_.ACES_solver
             info(1) = abs(eu(1)+eu(2));
             info(2) = 1.0e+8;
             %     return
@@ -370,7 +370,7 @@ end % end if useAIM and...
         dr.eigval = eig(G1pi);
         dr.rank=FL_RANK;
 
-        if options_.ACES_solver==1
+        if options_.ACES_solver
             betap=options_.planner_discount;
             sigma_cov=M_.Sigma_e;
             % get W - BY
@@ -440,7 +440,7 @@ end % end if useAIM and...
 
     catch
         lerror=lasterror;
-        if options_.ACES_solver==1
+        if options_.ACES_solver
             disp('Problem with using Part Info ACES solver:');
             error(lerror.message);
         else
@@ -451,4 +451,4 @@ end % end if useAIM and...
 
     % TODO:
     % if options_.loglinear == 1
-    % if exogenous deterministic variables
\ No newline at end of file
+    % if exogenous deterministic variables
diff --git a/matlab/print_bytecode_dynamic_model.m b/matlab/print_bytecode_dynamic_model.m
index af5eb011546ead3cb920b74b90e9fd979a32a876..09613f6e3582573e88eb9ade7297c5e8edf56f68 100644
--- a/matlab/print_bytecode_dynamic_model.m
+++ b/matlab/print_bytecode_dynamic_model.m
@@ -28,8 +28,8 @@ function print_bytecode_dynamic_model()
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 global options_
-if (options_.bytecode == 1)
+if options_.bytecode
     bytecode('print','dynamic');
 else
     disp('You have to use bytecode option in model command to use print_bytecode_dynamic_model');
-end
\ No newline at end of file
+end
diff --git a/matlab/print_bytecode_static_model.m b/matlab/print_bytecode_static_model.m
index ece12428b5237bd5ceb9a7293b2d727423b78a2a..5eb1427bd1f47de5e58288f0050774aa611be727 100644
--- a/matlab/print_bytecode_static_model.m
+++ b/matlab/print_bytecode_static_model.m
@@ -28,8 +28,8 @@ function print_bytecode_static_model()
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 global options_
-if (options_.bytecode == 1)
+if options_.bytecode
     bytecode('print','static');
 else
     disp('You have to use bytecode option in model command to use print_bytecode_static_model');
-end
\ No newline at end of file
+end
diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m
index 320fb272febfd39120672a735d63df25ae364cfb..7559b7c3cacbe57eeb81291486b964dba3d6ee21 100644
--- a/matlab/prior_posterior_statistics_core.m
+++ b/matlab/prior_posterior_statistics_core.m
@@ -226,9 +226,9 @@ for b=fpar:B
         %% Compute constant for observables
         if options_.prefilter == 1 %as mean is taken after log transformation, no distinction is needed here
             constant_part=repmat(mean_varobs',1,gend);
-        elseif options_.prefilter == 0 && options_.loglinear == 1 %logged steady state must be used
+        elseif options_.prefilter == 0 && options_.loglinear %logged steady state must be used
             constant_part=repmat(log(SteadyState(IdObs)),1,gend);
-        elseif options_.prefilter == 0 && options_.loglinear == 0 %unlogged steady state must be used
+        elseif options_.prefilter == 0 && ~options_.loglinear %unlogged steady state must be used
             constant_part=repmat(SteadyState(IdObs),1,gend);
         end
         %add trend to observables
diff --git a/matlab/ramsey_policy.m b/matlab/ramsey_policy.m
index bde0aaa1289e741a06f24dd8153d70749ba45a03..145aa56e6b61c53e81b67c91c343e7dc74d70db1 100644
--- a/matlab/ramsey_policy.m
+++ b/matlab/ramsey_policy.m
@@ -43,7 +43,7 @@ info = stoch_simul(var_list);
 
 oo_.steady_state = oo_.dr.ys;
 
-if options_.noprint == 0
+if ~options_.noprint
     disp_steady_state(M_,oo_)
     for i=M_.orig_endo_nbr:M_.endo_nbr
         if strmatch('mult_', M_.endo_names{i})
@@ -55,4 +55,4 @@ end
 
 oo_.planner_objective_value = evaluate_planner_objective(M_,options_,oo_);
 
-options_ = oldoptions;
\ No newline at end of file
+options_ = oldoptions;
diff --git a/matlab/set_default_plot_shock_decomposition_options.m b/matlab/set_default_plot_shock_decomposition_options.m
index 327c047cb9f04c6b94b3dc2f06201e7fc83d21f2..dc58373d774846af2c8a1a9fd06155977a1744c7 100644
--- a/matlab/set_default_plot_shock_decomposition_options.m
+++ b/matlab/set_default_plot_shock_decomposition_options.m
@@ -30,15 +30,15 @@ function options = set_default_plot_shock_decomposition_options(options)
 
 options.plot_shock_decomp.use_shock_groups = '';
 options.plot_shock_decomp.colormap = '';
-options.plot_shock_decomp.nodisplay = 0;
+options.plot_shock_decomp.nodisplay = false;
 options.plot_shock_decomp.graph_format = 'eps';
-options.plot_shock_decomp.detail_plot = 0;
-options.plot_shock_decomp.interactive = 0;
-options.plot_shock_decomp.screen_shocks = 0;
-options.plot_shock_decomp.steadystate = 0;
+options.plot_shock_decomp.detail_plot = false;
+options.plot_shock_decomp.interactive = false;
+options.plot_shock_decomp.screen_shocks = false;
+options.plot_shock_decomp.steadystate = false;
 options.plot_shock_decomp.type = '';
 options.plot_shock_decomp.fig_name = '';
-options.plot_shock_decomp.write_xls = 0;
+options.plot_shock_decomp.write_xls = false;
 options.plot_shock_decomp.realtime = 0; % 0 is standard; 1 is realtime
                                         % (pool/vintage); 2 is conditional
                                         % (pool/vintage); 3 is forecast
diff --git a/matlab/set_state_space.m b/matlab/set_state_space.m
index de8108f3273ad7898d44133777f2b1fb4ad94679..209a60b1f1010e003c8015e52d0c7baf4c98848e 100644
--- a/matlab/set_state_space.m
+++ b/matlab/set_state_space.m
@@ -69,7 +69,7 @@ else
     both_var = [];
     stat_var = setdiff([1:endo_nbr]',fwrd_var);
 end
-if DynareOptions.block == 1
+if DynareOptions.block
     order_var = DynareModel.block_structure.variable_reordered;
 else
     order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)];
diff --git a/matlab/solve1.m b/matlab/solve1.m
index 5ea2651bfd869f154486f4edcf0dfcc3f0d91123..ba3c263b95817bf98fae62e59b8d4dc5bb33dfcb 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -6,8 +6,8 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,gstep,tolf,tolx,maxit,deb
 %    x:               guess values
 %    j1:              equations index for which the model is solved
 %    j2:              unknown variables index
-%    jacobian_flag=1: jacobian given by the 'func' function
-%    jacobian_flag=0: jacobian obtained numerically
+%    jacobian_flag=true: jacobian given by the 'func' function
+%    jacobian_flag=false: jacobian obtained numerically
 %    gstep            increment multiplier in numercial derivative
 %                     computation
 %    tolf             tolerance for residuals
diff --git a/matlab/steady.m b/matlab/steady.m
index ace9a5771dcf3ecaaff5c4a277c0ff76223015aa..ba8beca7df5cf4a322c42d9a7d612992759729c2 100644
--- a/matlab/steady.m
+++ b/matlab/steady.m
@@ -80,11 +80,11 @@ end
 [oo_.steady_state,M_.params,info] = steady_(M_,options_,oo_);
 
 if info(1) == 0
-    if options_.noprint == 0
+    if ~options_.noprint
         disp_steady_state(M_,oo_);
     end
 else
-    if options_.noprint == 0
+    if ~options_.noprint
         if ~isempty(oo_.steady_state)
             resid;
         else
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 3be3c827648f7cffcfd0a50b8053a082d749ed2d..4e7d6a41a8f24cd441949ad2e3c8c2bdcabb1329 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -45,7 +45,7 @@ if isempty(options_.qz_criterium)
     options_.qz_criterium = 1+1e-6;
 end
 
-if options_.partial_information == 1 || options_.ACES_solver == 1
+if options_.partial_information || options_.ACES_solver
     PI_PCL_solver = 1;
     if options_.order ~= 1
         warning('stoch_simul:: forcing order=1 since you are using partial_information or ACES solver')
@@ -182,7 +182,7 @@ if options_.periods > 0 && ~PI_PCL_solver
     end
 end
 
-if options_.nomoments == 0
+if ~options_.nomoments
     if PI_PCL_solver
         PCL_Part_info_moments(0, PCL_varobs, oo_.dr, i_var);
     elseif options_.periods == 0
@@ -266,7 +266,7 @@ if options_.irf
                     end
                 end
             end
-            if options_.nograph == 0
+            if ~options_.nograph
                 number_of_plots_to_draw = size(irfs,1);
                 [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw);
                 if nbplt == 0
@@ -379,7 +379,7 @@ if options_.irf
     end
 end
 
-if options_.SpectralDensity.trigger == 1
+if options_.SpectralDensity.trigger
     [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list);
 end
 
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index 9f08b16010a4fd0435a064c7a97ce4f8c3741464..2c2591f9689a9ad0a0263d48b08541da25a813a5 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -56,7 +56,7 @@ if options_.order>2 && ~options_.k_order_solver
     error('You need to set k_order_solver for order>2')
 end
 
-if (options_.aim_solver == 1) && (local_order > 1)
+if options_.aim_solver && (local_order > 1)
     error('Option "aim_solver" is incompatible with order >= 2')
 end
 
@@ -255,7 +255,7 @@ elseif options_.risky_steadystate
     options_.order = orig_order;
 else
     % If required, use AIM solver if not check only
-    if (options_.aim_solver == 1) && (task == 0)
+    if options_.aim_solver && (task == 0)
         [dr,info] = AIM_first_order_solver(jacobia_,M_,dr,options_.qz_criterium);
 
     else  % use original Dynare solver
diff --git a/matlab/store_smoother_results.m b/matlab/store_smoother_results.m
index bad4fe723d2c4faeb6454de47d12291d576842f6..cc529fab24c183b9a3d991c8fb92b62459d28bf6 100644
--- a/matlab/store_smoother_results.m
+++ b/matlab/store_smoother_results.m
@@ -102,9 +102,9 @@ end
 %% Compute constant for observables
 if options_.prefilter == 1 %as mean is taken after log transformation, no distinction is needed here
     constant_part=repmat(dataset_info.descriptive.mean',1,gend);
-elseif options_.prefilter == 0 && options_.loglinear == 1 %logged steady state must be used
+elseif options_.prefilter == 0 && options_.loglinear %logged steady state must be used
     constant_part=repmat(log(ys(bayestopt_.mfys)),1,gend);
-elseif options_.prefilter == 0 && options_.loglinear == 0 %unlogged steady state must be used
+elseif options_.prefilter == 0 && ~options_.loglinear %unlogged steady state must be used
     constant_part=repmat(ys(bayestopt_.mfys),1,gend);
 end
 
@@ -134,9 +134,9 @@ i_endo_in_dr_matrices=bayestopt_.smoother_var_list(i_endo_in_bayestopt_smoother_
 if ~isempty(options_.nk) && options_.nk ~= 0
     %% Compute constant
     i_endo_declaration_order = oo_.dr.order_var(i_endo_in_dr_matrices); %get indices of smoothed variables in name vector
-    if  options_.loglinear == 1 %logged steady state must be used
+    if  options_.loglinear %logged steady state must be used
         constant_all_variables=repmat(log(ys(i_endo_declaration_order))',[length(options_.filter_step_ahead),1,gend+max(options_.filter_step_ahead)]);
-    elseif options_.loglinear == 0 %unlogged steady state must be used
+    elseif ~options_.loglinear %unlogged steady state must be used
         constant_all_variables=repmat((ys(i_endo_declaration_order))',[length(options_.filter_step_ahead),1,gend+max(options_.filter_step_ahead)]);
     end
     % add constant
diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m
index 3cde22f314cc760519ad40a94fcf1c77261e4fc4..4b82390437feb0b5fda63a27af0f7d76dacb022d 100644
--- a/matlab/th_autocovariances.m
+++ b/matlab/th_autocovariances.m
@@ -92,7 +92,7 @@ nspred = M_.nspred;
 nstatic = M_.nstatic;
 
 nx = size(ghx,2);
-if options_.block == 0
+if ~options_.block
     %order_var = dr.order_var;
     inv_order_var = dr.inv_order_var;
     kstate = dr.kstate;
@@ -123,7 +123,7 @@ end
 b = ghu1*M_.Sigma_e*ghu1';
 
 
-if options_.block == 0
+if ~options_.block
     ipred = nstatic+(1:nspred)';
 else
     ipred = dr.state_var;
@@ -135,7 +135,7 @@ end
 % HP filtering, this mean correction is computed *before* filtering)
 if local_order == 2 || options_.hp_filter == 0
     [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
-    if options_.block == 0
+    if ~options_.block
         iky = inv_order_var(ivar);
     else
         iky = ivar;
diff --git a/matlab/trust_region.m b/matlab/trust_region.m
index 13f11bd3da1cab3e6947bce3400030f359e10486..f4b445cc7be2d91e4e6bfa972c9baf2c4bee894d 100644
--- a/matlab/trust_region.m
+++ b/matlab/trust_region.m
@@ -7,8 +7,8 @@ function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tol
 %    x0:              guess values
 %    j1:              equations index for which the model is solved
 %    j2:              unknown variables index
-%    jacobian_flag=1: jacobian given by the 'func' function
-%    jacobian_flag=0: jacobian obtained numerically
+%    jacobian_flag=true: jacobian given by the 'func' function
+%    jacobian_flag=false: jacobian obtained numerically
 %    gstep            increment multiplier in numercial derivative
 %                     computation
 %    tolf             tolerance for residuals
diff --git a/preprocessor b/preprocessor
index d9f7ac4c9be88f1e6d8b953c0c71195e1fa31f8c..78583135dfbc0e87a2bca83481d2d7939c3467e9 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit d9f7ac4c9be88f1e6d8b953c0c71195e1fa31f8c
+Subproject commit 78583135dfbc0e87a2bca83481d2d7939c3467e9