diff --git a/matlab/display_problematic_vars_Jacobian.m b/matlab/display_problematic_vars_Jacobian.m
index 0cb42b581..b77c60804 100644
--- a/matlab/display_problematic_vars_Jacobian.m
+++ b/matlab/display_problematic_vars_Jacobian.m
@@ -36,6 +36,11 @@ function []=display_problematic_vars_Jacobian(problemrow,problemcol,M_,x,type,ca
 if nargin<6
     caller_string='';
 end
+aux_vars_type = [M_.aux_vars.type];
+aux_eq_nbr = max(find(aux_vars_type == 6)); 
+if isempty(aux_eq_nbr)
+    aux_eq_nbr=0;
+end
 if strcmp(type,'dynamic')
     for ii=1:length(problemrow)
         [var_row,var_index]=find(M_.lead_lag_incidence==problemcol(ii));
@@ -47,10 +52,33 @@ if strcmp(type,'dynamic')
             type_string='lead of';
         end
         if var_index<=M_.orig_endo_nbr
-            fprintf('Derivative of Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',problemrow(ii),type_string,deblank(M_.endo_names(var_index,:)),deblank(M_.endo_names(var_index,:)),x(var_index))
+            if problemrow(ii)<=aux_eq_nbr
+                eq_nbr=problemrow(ii);
+                fprintf('Derivative of Auxiliary Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',eq_nbr,type_string,deblank(M_.endo_names(var_index,:)),deblank(M_.endo_names(var_index,:)),x(var_index))
+            else
+                eq_nbr=problemrow(ii)-aux_eq_nbr;
+                fprintf('Derivative of Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',eq_nbr,type_string,deblank(M_.endo_names(var_index,:)),deblank(M_.endo_names(var_index,:)),x(var_index))
+            end
         else %auxiliary vars
-            orig_var_index=M_.aux_vars(1,var_index-M_.orig_endo_nbr).orig_index;
-            fprintf('Derivative of Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',problemrow(ii),type_string,deblank(M_.endo_names(orig_var_index,:)),deblank(M_.endo_names(orig_var_index,:)),x(orig_var_index))            
+            if M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).type ==6 %Ramsey Lagrange Multiplier 
+                if problemrow(ii)<=aux_eq_nbr
+                    eq_nbr=problemrow(ii);
+                    fprintf('Derivative of Auxiliary Equation %d with respect to %s of Langrange multiplier of equation %s (initial value: %g) \n',eq_nbr,type_string,M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).eq_nbr,x(problemcol(ii)))            
+                else
+                    eq_nbr=problemrow(ii)-aux_eq_nbr;
+                    fprintf('Derivative of Equation %d with respect to %s of Langrange multiplier of equation %s (initial value: %g) \n',eq_nbr,type_string,M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).eq_nbr,x(problemcol(ii)))            
+                end
+            else
+                if problemrow(ii)<=aux_eq_nbr
+                    eq_nbr=problemrow(ii);
+                    orig_var_index=M_.aux_vars(1,var_index-M_.orig_endo_nbr).orig_index;
+                    fprintf('Derivative of Auxiliary Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',eq_nbr,type_string,deblank(M_.endo_names(orig_var_index,:)),deblank(M_.endo_names(orig_var_index,:)),x(orig_var_index))            
+                else
+                    eq_nbr=problemrow(ii)-aux_eq_nbr;
+                    orig_var_index=M_.aux_vars(1,var_index-M_.orig_endo_nbr).orig_index;
+                    fprintf('Derivative of Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',eq_nbr,type_string,deblank(M_.endo_names(orig_var_index,:)),deblank(M_.endo_names(orig_var_index,:)),x(orig_var_index))            
+                end
+            end
         end    
     end
     fprintf('\n%s  The problem most often occurs, because a variable with\n',caller_string)
@@ -59,10 +87,33 @@ if strcmp(type,'dynamic')
 elseif strcmp(type,'static')
     for ii=1:length(problemrow)
         if problemcol(ii)<=M_.orig_endo_nbr
-            fprintf('Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',problemrow(ii),deblank(M_.endo_names(problemcol(ii),:)),deblank(M_.endo_names(problemcol(ii),:)),x(problemcol(ii)))
+            if problemrow(ii)<=aux_eq_nbr
+                eq_nbr=problemrow(ii);
+                fprintf('Derivative of Auxiliary Equation %d with respect to Variable %s  (initial value of %s: %g) \n',eq_nbr,deblank(M_.endo_names(problemcol(ii),:)),deblank(M_.endo_names(problemcol(ii),:)),x(problemcol(ii)))
+            else
+                eq_nbr=problemrow(ii)-aux_eq_nbr;
+                fprintf('Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',eq_nbr,deblank(M_.endo_names(problemcol(ii),:)),deblank(M_.endo_names(problemcol(ii),:)),x(problemcol(ii)))
+            end
         else %auxiliary vars
-            orig_var_index=M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).orig_index;
-            fprintf('Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',problemrow(ii),deblank(M_.endo_names(orig_var_index,:)),deblank(M_.endo_names(orig_var_index,:)),x(problemcol(ii)))            
+            if M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).type ==6 %Ramsey Lagrange Multiplier 
+                if problemrow(ii)<=aux_eq_nbr
+                    eq_nbr=problemrow(ii);
+                    fprintf('Derivative of Auxiliary Equation %s with respect to Lagrange multiplier of equation %d (initial value: %g) \n',eq_nbr,M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).eq_nbr,x(problemcol(ii)))            
+                else
+                    eq_nbr=problemrow(ii)-aux_eq_nbr;
+                    fprintf('Derivative of Equation %s with respect to Lagrange multiplier of equation %d (initial value: %g) \n',eq_nbr,M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).eq_nbr,x(problemcol(ii)))            
+                end
+            else
+            if problemrow(ii)<=aux_eq_nbr
+                eq_nbr=problemrow(ii);
+                orig_var_index=M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).orig_index;
+                fprintf('Derivative of Auxiliary Equation %d with respect to Variable %s  (initial value of %s: %g) \n',eq_nbr,deblank(M_.endo_names(orig_var_index,:)),deblank(M_.endo_names(orig_var_index,:)),x(problemcol(ii)))            
+            else
+                eq_nbr=problemrow(ii)-aux_eq_nbr;
+                orig_var_index=M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).orig_index;
+                fprintf('Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',eq_nbr,deblank(M_.endo_names(orig_var_index,:)),deblank(M_.endo_names(orig_var_index,:)),x(problemcol(ii)))            
+            end
+            end
         end
     end
     fprintf('\n%s  The problem most often occurs, because a variable with\n',caller_string)
diff --git a/matlab/dyn_ramsey_static.m b/matlab/dyn_ramsey_static.m
index 32eaa745b..96463aeb5 100644
--- a/matlab/dyn_ramsey_static.m
+++ b/matlab/dyn_ramsey_static.m
@@ -14,7 +14,7 @@ function [steady_state,params,check] = dyn_ramsey_static(x,M,options_,oo)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2003-2012 Dynare Team
+% Copyright (C) 2003-2014 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -32,20 +32,18 @@ function [steady_state,params,check] = dyn_ramsey_static(x,M,options_,oo)
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 
-steady_state = [];
 params = M.params;
 check = 0;
-options_.steadystate.nocheck = 1;
+options_.steadystate.nocheck = 1; %disable checking because Lagrange multipliers are not accounted for in evaluate_steady_state_file
 % dyn_ramsey_static_1 is a subfunction
 nl_func = @(x) dyn_ramsey_static_1(x,M,options_,oo);
 
 % check_static_model is a subfunction
-if check_static_model(oo.steady_state,M,options_,oo)
+if check_static_model(oo.steady_state,M,options_,oo) && ~options_.steadystate_flag
     steady_state = oo.steady_state;
     return
 elseif options_.steadystate_flag
     k_inst = [];
-    instruments = options_.instruments;
     inst_nbr = size(options_.instruments,1);
     for i = 1:inst_nbr
         k_inst = [k_inst; strmatch(options_.instruments(i,:), ...
@@ -54,16 +52,16 @@ elseif options_.steadystate_flag
     if inst_nbr == 1
         inst_val = csolve(nl_func,oo.steady_state(k_inst),'',options_.solve_tolf,100);
     else
-        [inst_val,info1] = dynare_solve(nl_func,ys(k_inst),0);
+        [inst_val,check] = dynare_solve(nl_func,ys(k_inst),0);
     end
     ys(k_inst) = inst_val;
     exo_ss = [oo.exo_steady_state oo.exo_det_steady_state];
     [xx,params,check] = evaluate_steady_state_file(ys,exo_ss,M,options_);
-    [junk,jun,steady_state] = nl_func(inst_val);
+    [junk,junk,steady_state] = nl_func(inst_val);
 else
     n_var = M.orig_endo_nbr;
     xx = oo.steady_state(1:n_var);
-    [xx,info1] = dynare_solve(nl_func,xx,0);
+    [xx,check] = dynare_solve(nl_func,xx,0);
     [junk,junk,steady_state] = nl_func(xx);
 end
 
@@ -78,24 +76,14 @@ mult = [];
 params = M.params;
 endo_nbr = M.endo_nbr;
 endo_names = M.endo_names;
-exo_nbr = M.exo_nbr;
 orig_endo_nbr = M.orig_endo_nbr;
 aux_vars_type = [M.aux_vars.type];
 orig_endo_aux_nbr = orig_endo_nbr + min(find(aux_vars_type == 6)) - 1; 
 orig_eq_nbr = M.orig_eq_nbr;
 inst_nbr = orig_endo_aux_nbr - orig_eq_nbr;
 % indices of Lagrange multipliers
-i_mult = [orig_endo_aux_nbr+(1:orig_eq_nbr)]';
 fname = M.fname;
-max_lead = M.maximum_lead;
-max_lag = M.maximum_lag;
 
-% indices of all endogenous variables
-i_endo = [1:endo_nbr]';
-% indices of endogenous variable except instruments
-% i_inst = M.instruments;
-% lead_lag incidence matrix
-i_lag = M.lead_lag_incidence;
 
 if options_.steadystate_flag
     k_inst = [];
@@ -112,7 +100,7 @@ if options_.steadystate_flag
 end
 
 xx = zeros(endo_nbr,1);
-xx(1:length(x)) = x;
+xx(1:M.orig_endo_nbr) = x(1:M.orig_endo_nbr); %take care of steady state file that will also return multipliers
 % setting steady state of auxiliary variables
 % that depends on original endogenous variables
 if any([M.aux_vars.type] ~= 6)
@@ -130,7 +118,6 @@ end
 % value and Jacobian of objective function
 ex = zeros(1,M.exo_nbr);
 [U,Uy,Uyy] = feval([fname '_objective_static'],x,ex, params);
-Uy = Uy';
 Uyy = reshape(Uyy,endo_nbr,endo_nbr);
 
 % set multipliers and auxiliary variables that
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index dcaa56133..53d842de7 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -52,7 +52,119 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
     end
 
     if options.ramsey_policy
-        [ys,params] = dyn_ramsey_static(ys_init,M,options,oo);
+        if steadystate_flag
+            % explicit steady state file
+            [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M, ...
+                                                           options);
+            %test whether it solves model conditional on the instruments
+            [resids] = evaluate_static_model(ys,exo_ss,params,M,options); 
+            n_multipliers=M.endo_nbr-M.orig_endo_nbr;
+            nan_indices=find(isnan(resids(n_multipliers+1:end)));
+
+            if ~isempty(nan_indices)
+                fprintf('\nevaluate_steady_state: The steady state file computation for the Ramsey problem resulted in NaNs.\n')
+                fprintf('evaluate_steady_state: The steady state was computed conditional on the following initial instrument values: \n') 
+                for ii = 1:size(options.instruments,1);
+                    fprintf('\t %s \t %f \n',options.instruments(ii,:),ys_init(strmatch(options.instruments(ii,:),M.endo_names,'exact')))
+                end
+                fprintf('evaluate_steady_state: The problem occured in the following equations: \n')
+                fprintf('\t Equation(s): ')
+                for ii=1:length(nan_indices)
+                        fprintf('%d, ',nan_indices(ii));
+                end
+                skipline();
+                fprintf('evaluate_steady_state: If those initial values are not admissable, change them using an initval-block.\n') 
+                skipline(2);
+                check=1;
+                info(1) = 84;
+                info(2) = resids'*resids;           
+                return;
+            end
+            if max(abs(resids(n_multipliers+1:end))) > options.dynatol.f %does it solve for all variables except for the Lagrange multipliers
+                fprintf('\nevaluate_steady_state: The steady state file does not solve the steady state for the Ramsey problem.\n')
+                fprintf('evaluate_steady_state: Conditional on the following instrument values: \n') 
+                for ii = 1:size(options.instruments,1);
+                    fprintf('\t %s \t %f \n',options.instruments(ii,:),ys_init(strmatch(options.instruments(ii,:),M.endo_names,'exact')))
+                end
+                fprintf('evaluate_steady_state: the following equations have non-zero residuals: \n')                
+                for ii=n_multipliers+1:M.endo_nbr
+                    if abs(resids(ii)) > options.dynatol.f/100
+                        fprintf('\t Equation number %d: %f\n',ii-n_multipliers, resids(ii))
+                    end
+                end
+                skipline(2);
+                info(1) = 85;
+                info(2) = resids'*resids;           
+                return;
+            end
+        end
+        %either if no steady state file or steady state file without problems
+        [ys,params,info] = dyn_ramsey_static(ys_init,M,options,oo);
+        if info
+           info=81;%case should not happen
+           return;
+        end
+        %check whether steady state really solves the model 
+        [resids] = evaluate_static_model(ys,exo_ss,params,M,options); 
+
+        n_multipliers=M.endo_nbr-M.orig_endo_nbr;
+        nan_indices_multiplier=find(isnan(resids(1:n_multipliers)));
+        nan_indices=find(isnan(resids(n_multipliers+1:end)));
+
+        if ~isempty(nan_indices)
+            fprintf('\nevaluate_steady_state: The steady state computation for the Ramsey problem resulted in NaNs.\n')
+            fprintf('evaluate_steady_state: The steady state computation resulted in the following instrument values: \n') 
+            for i = 1:size(options.instruments,1);
+                fprintf('\t %s \t %f \n',options.instruments(i,:),ys(strmatch(options.instruments(i,:),M.endo_names,'exact')))
+            end
+            fprintf('evaluate_steady_state: The problem occured in the following equations: \n')
+            fprintf('\t Equation(s): ')
+            for ii=1:length(nan_indices)
+                    fprintf('%d, ',nan_indices(ii));
+            end
+            skipline();
+            info(1) = 82;
+            return;
+        end
+
+        if ~isempty(nan_indices_multiplier)
+            fprintf('\nevaluate_steady_state: The steady state computation for the Ramsey problem resulted in NaNs in the auxiliary equations.\n')
+            fprintf('evaluate_steady_state: The steady state computation resulted in the following instrument values: \n') 
+            for i = 1:size(options.instruments,1);
+                fprintf('\t %s \t %f \n',options.instruments(i,:),ys(strmatch(options.instruments(i,:),M.endo_names,'exact')))
+            end
+            fprintf('evaluate_steady_state: The problem occured in the following equations: \n')
+            fprintf('\t Auxiliary equation(s): ')
+            for ii=1:length(nan_indices_multiplier)
+                    fprintf('%d, ',nan_indices_multiplier(ii));
+            end
+            skipline();
+            info(1) = 83;
+            return;
+        end
+
+        if max(abs(resids)) > options.dynatol.f %does it solve for all variables including the auxiliary ones
+            fprintf('\nevaluate_steady_state: The steady state for the Ramsey problem could not be computed.\n')
+            fprintf('evaluate_steady_state: The steady state computation stopped with the following instrument values:: \n') 
+            for i = 1:size(options.instruments,1);
+                fprintf('\t %s \t %f \n',options.instruments(i,:),ys_init(strmatch(options.instruments(i,:),M.endo_names,'exact')))
+            end
+            fprintf('evaluate_steady_state: The following equations have non-zero residuals: \n')                
+            for ii=1:n_multipliers
+                if abs(resids(ii)) > options.dynatol.f/100
+                    fprintf('\t Auxiliary Ramsey equation number %d: %f\n',ii, resids(ii))
+                end
+            end
+            for ii=n_multipliers+1:M.endo_nbr
+                if abs(resids(ii)) > options.dynatol.f/100
+                    fprintf('\t Equation number %d: %f\n',ii-n_multipliers, resids(ii))
+                end
+            end
+            skipline(2);
+            info(1) = 81;
+            info(2) = resids'*resids;           
+            return;
+        end       
     elseif steadystate_flag
         % explicit steady state file
         [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M, ...
@@ -109,13 +221,8 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
     end
 
     if check
-        if options.steadystate_flag
-            info(1)= 19;
-            resid = check1 ;
-        else
-            info(1)= 20;
-            resid = evaluate_static_model(ys_init,exo_ss,params,M,options);
-        end
+        info(1)= 20;
+        resid = evaluate_static_model(ys_init,exo_ss,params,M,options);
         info(2) = resid'*resid ;
         if isnan(info(2))
             info(1)=22;
@@ -161,3 +268,5 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
         return
     end
 
+    
+
diff --git a/matlab/print_info.m b/matlab/print_info.m
index 2aded2e6b..f3d0c6ac2 100644
--- a/matlab/print_info.m
+++ b/matlab/print_info.m
@@ -120,6 +120,18 @@ if ~noprint
       case 72
         error(['Calibrated covariance of the measurement errors implies correlation larger than  +-1.']);
         % Aim Code Conversions by convertAimCodeToInfo.m
+      case 81
+        error(['Ramsey: The solution to the static first order conditions for optimal policy could not be found. Either the model' ...
+               ' doesn''t have a steady state, there are an infinity of steady states, ' ...
+               ' or the guess values are too far from the solution']);
+      case 82
+        error(['Ramsey: The steady state computation resulted in NaN in the static first order conditions for optimal policy']);
+      case 83
+        error(['Ramsey: The steady state computation resulted in NaN in the auxiliary equations for optimal policy']);
+      case 84
+        error(['Ramsey: The steady state file computation for the Ramsey problem resulted in NaNs at the initial values of the instruments']);
+      case 85
+        error(['Ramsey: The steady state file does not solve the static first order conditions conditional on the instruments.']);
       case 102
         error('Aim: roots not correctly computed by real_schur');
       case 103
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index ca2051329..0104c3682 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -1523,7 +1523,8 @@ DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
                     << "%" << endl
                     << "% Outputs:" << endl
                     << "%   residual  [M_.endo_nbr by 1] double    vector of residuals of the dynamic model equations in order of " << endl
-                    << "%                                          declaration of the equations" << endl
+                    << "%                                          declaration of the equations." << endl
+					<< "%                                          Dynare may prepend auxiliary equations, see M_.aux_vars" << endl
                     << "%   g1        [M_.endo_nbr by #dynamic variables] double    Jacobian matrix of the dynamic model equations;" << endl
                     << "%                                                           rows: equations in order of declaration" << endl
                     << "%                                                           columns: variables in order stored in M_.lead_lag_incidence" << endl
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index 18be48f11..09fb9261b 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -1159,7 +1159,8 @@ StaticModel::writeStaticMFile(const string &func_name) const
          << "%" << endl
          << "% Outputs:" << endl
          << "%   residual  [M_.endo_nbr by 1] double    vector of residuals of the static model equations " << endl
-         << "%                                          in order of declaration of the equations" << endl
+         << "%                                          in order of declaration of the equations." << endl
+         << "%                                          Dynare may prepend auxiliary equations, see M_.aux_vars" << endl
          << "%   g1        [M_.endo_nbr by M_.endo_nbr] double    Jacobian matrix of the static model equations;" << endl
          << "%                                                       columns: variables in declaration order" << endl
          << "%                                                       rows: equations in order of declaration" << endl
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 257af1500..0147e1621 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -26,6 +26,9 @@ MODFILES = \
 	optimal_policy/nk_ramsey_expectation.mod \
 	optimal_policy/nk_ramsey_expectation_a.mod \
 	optimal_policy/mult_elimination_test.mod \
+	optimal_policy/Ramsey/ramsey_ex_initval.mod \
+	optimal_policy/Ramsey/ramsey_ex.mod \
+	optimal_policy/Ramsey/ramsey_ex_initval_AR2.mod \
 	discretionary_policy/dennis_1.mod \
 	ramst_initval_file.mod \
 	ramst_normcdf_and_friends.mod \
@@ -184,7 +187,8 @@ MODFILES = \
 	reporting/example1.mod
 
 XFAIL_MODFILES = ramst_xfail.mod \
-	estim_param_in_shock_value.mod
+	estim_param_in_shock_value.mod \
+	optimal_policy/Ramsey/ramsey_ex_wrong_ss_file.mod \
 
 # Dependencies
 example1_use_dll.m.trs: example1.m.trs
@@ -346,7 +350,9 @@ EXTRA_DIST = \
 	filter_step_ahead/fs2000_filter_step_ahead_bayesian_steadystate.m \
 	filter_step_ahead/fs2000_filter_step_ahead_ML_steadystate.m \
 	loglinear/results_exp.mat \
-	smoother2histval/fsdat_simul.m
+	smoother2histval/fsdat_simul.m \
+	optimal_policy/Ramsey/find_c.m \
+	optimal_policy/Ramsey/oo_ramsey_policy_initval.mat \
 
 TARGETS =
 
diff --git a/tests/optimal_policy/Ramsey/find_c.m b/tests/optimal_policy/Ramsey/find_c.m
new file mode 100644
index 000000000..badafbf98
--- /dev/null
+++ b/tests/optimal_policy/Ramsey/find_c.m
@@ -0,0 +1,8 @@
+function c = find_c(c0,pai,beta,epsilon,phi,gamma,omega)
+c = csolve(@nk_ss,c0,[],1e-8,100,pai,beta,epsilon,phi,gamma,omega);
+end
+
+function r = nk_ss(c,pai,beta,epsilon,phi,gamma,omega)
+r = pai*(pai-1)/c - beta*pai*(pai-1)/c-epsilon*phi*(c+(omega/2)*(pai-1)^2)^(gamma+1)/omega +(c+(omega/2)*(pai-1)^2)*(epsilon-1)/(omega*c);
+
+end
\ No newline at end of file
diff --git a/tests/optimal_policy/Ramsey/oo_ramsey_policy_initval.mat b/tests/optimal_policy/Ramsey/oo_ramsey_policy_initval.mat
new file mode 100644
index 000000000..ff829dd9d
Binary files /dev/null and b/tests/optimal_policy/Ramsey/oo_ramsey_policy_initval.mat differ
diff --git a/tests/optimal_policy/Ramsey/ramsey_ex.mod b/tests/optimal_policy/Ramsey/ramsey_ex.mod
new file mode 100644
index 000000000..d28882f97
--- /dev/null
+++ b/tests/optimal_policy/Ramsey/ramsey_ex.mod
@@ -0,0 +1,63 @@
+/* Mod file tests the correctness of the Ramsey command when used together with a steady state file by
+ * - checking whether the results coincide with the ones when used with an initval block
+ * - checking whether between stoch_simul and ramsey_planner are consistent
+ * The example is taken from Juillard, Michel (2011): User manual for optimal policy package, 
+ * MONFISPOL FP7 project SSH-225149, Deliverable 1.1.2
+*/
+
+var pai, c, n, r, a;
+varexo u;
+parameters beta, rho, epsilon, omega, phi, gamma;
+
+beta=0.99;
+gamma=3;
+omega=17;
+epsilon=8;
+phi=1;
+rho=0.95;
+
+model;
+a = rho*a(-1)+u;
+1/c = beta*r/(c(+1)*pai(+1));
+pai*(pai-1)/c = beta*pai(+1)*(pai(+1)-1)/c(+1)+epsilon*phi*n^(gamma+1)/omega -exp(a)*n*(epsilon-1)/(omega*c);
+exp(a)*n = c+(omega/2)*(pai-1)^2;
+end;
+
+initval;
+r=1;
+end;
+
+steady_state_model;
+a = 0;
+pai = beta*r;
+c = find_c(0.96,pai,beta,epsilon,phi,gamma,omega);
+n = c+(omega/2)*(pai-1)^2;
+end;
+
+shocks;
+var u; stderr 0.008;
+end;
+
+planner_objective(ln(c)-phi*((n^(1+gamma))/(1+gamma)));
+ramsey_policy(planner_discount=0.99,order=1,instruments=(r));
+oo_ramsey_policy_steady_state_file=oo_;
+stoch_simul(periods=0, order=1, irf=25, nograph);
+if max(abs((oo_ramsey_policy_steady_state_file.steady_state-oo_.steady_state)))>1e-5 ...
+    || max(abs(oo_ramsey_policy_steady_state_file.dr.ys-oo_.dr.ys))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_steady_state_file.dr.ghx-oo_.dr.ghx)))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_steady_state_file.dr.ghu-oo_.dr.ghu)))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_steady_state_file.dr.Gy-oo_.dr.Gy)))>1e-5 ...
+    || max(abs((oo_ramsey_policy_steady_state_file.planner_objective_value-oo_.planner_objective_value)))>1e-5
+    error('Running stoch_simul after ramsey_policy leads to inconsistent results')
+end
+
+load oo_ramsey_policy_initval;
+
+if max(abs((oo_ramsey_policy_initval.steady_state-oo_.steady_state)))>1e-5 ...
+    || max(abs(oo_ramsey_policy_initval.dr.ys-oo_.dr.ys))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_initval.dr.ghx-oo_.dr.ghx)))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_initval.dr.ghu-oo_.dr.ghu)))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_initval.dr.Gy-oo_.dr.Gy)))>1e-5 ...
+    || max(abs((oo_ramsey_policy_initval.planner_objective_value-oo_.planner_objective_value)))>1e-5
+    error('Initval and steady state file yield different results')
+end
diff --git a/tests/optimal_policy/Ramsey/ramsey_ex_initval.mod b/tests/optimal_policy/Ramsey/ramsey_ex_initval.mod
new file mode 100644
index 000000000..170b956f9
--- /dev/null
+++ b/tests/optimal_policy/Ramsey/ramsey_ex_initval.mod
@@ -0,0 +1,52 @@
+/* Mod file tests the functioning of the Ramsey command when used together with an initval-block
+ * - Tests whether subsequent calls to ramsey_policy and stoch_simul are possible
+ * The example is taken from Juillard, Michel (2011): User manual for optimal policy package, 
+ * MONFISPOL FP7 project SSH-225149, Deliverable 1.1.2
+*/
+
+var pai, c, n, r, a;
+varexo u;
+parameters beta, rho, epsilon, omega, phi, gamma;
+
+beta=0.99;
+gamma=3;
+omega=17;
+epsilon=8;
+phi=1;
+rho=0.95;
+
+model;
+a = rho*a(-1)+u;
+1/c = beta*r/(c(+1)*pai(+1));
+pai*(pai-1)/c = beta*pai(+1)*(pai(+1)-1)/c(+1)+epsilon*phi*n^(gamma+1)/omega -exp(a)*n*(epsilon-1)/(omega*c);
+exp(a)*n = c+(omega/2)*(pai-1)^2;
+end;
+
+initval;
+r=1;
+end;
+
+initval;
+a = 0;
+pai = beta;
+c = 0.9665;
+n = 0.9673;
+end;
+
+shocks;
+var u; stderr 0.008;
+end;
+
+planner_objective(ln(c)-phi*((n^(1+gamma))/(1+gamma)));
+ramsey_policy(planner_discount=0.99,order=1,instruments=(r));
+oo_ramsey_policy_initval=oo_;
+%save oo_ramsey_policy_initval oo_ramsey_policy_initval
+stoch_simul(periods=0, order=1, irf=25, nograph);
+if max(abs((oo_ramsey_policy_initval.steady_state-oo_.steady_state)))>1e-5 ...
+    || max(abs(oo_ramsey_policy_initval.dr.ys-oo_.dr.ys))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_initval.dr.ghx-oo_.dr.ghx)))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_initval.dr.ghu-oo_.dr.ghu)))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_initval.dr.Gy-oo_.dr.Gy)))>1e-5 ...
+    || max(abs((oo_ramsey_policy_initval.planner_objective_value-oo_.planner_objective_value)))>1e-5
+    error('Running stoch_simul after ramsey_policy leads to inconsistent results')
+end
diff --git a/tests/optimal_policy/Ramsey/ramsey_ex_initval_AR2.mod b/tests/optimal_policy/Ramsey/ramsey_ex_initval_AR2.mod
new file mode 100644
index 000000000..c10204922
--- /dev/null
+++ b/tests/optimal_policy/Ramsey/ramsey_ex_initval_AR2.mod
@@ -0,0 +1,50 @@
+/* Mod file tests the correctness of the Ramsey command when Auxiliary variables for AR2 are used
+ * The example is adapted from Juillard, Michel (2011): User manual for optimal policy package, 
+ * MONFISPOL FP7 project SSH-225149, Deliverable 1.1.2
+*/
+
+var pai, c, n, r, a;
+varexo u;
+parameters beta, rho, epsilon, omega, phi, gamma;
+
+beta=0.99;
+gamma=3;
+omega=17;
+epsilon=8;
+phi=1;
+rho=0.95;
+
+model;
+log(a) = rho*log(a(-2))+u;
+1/c = beta*r/(c(+1)*pai(+1));
+pai*(pai-1)/c = beta*pai(+1)*(pai(+1)-1)/c(+1)+epsilon*phi*n^(gamma+1)/omega -a*n*(epsilon-1)/(omega*c);
+a*n = c+(omega/2)*(pai-1)^2;
+end;
+
+initval;
+r=1;
+end;
+
+initval;
+a = 1;
+pai = beta;
+c = 0.9665;
+n = 0.9673;
+end;
+
+shocks;
+var u; stderr 0.008;
+end;
+
+planner_objective(ln(c)-phi*((n^(1+gamma))/(1+gamma)));
+ramsey_policy(planner_discount=0.99,order=1,instruments=(r));
+oo_ramsey_policy_initval_AR2=oo_;
+stoch_simul(periods=0, order=1, irf=25, nograph);
+if max(abs((oo_ramsey_policy_initval_AR2.steady_state-oo_.steady_state)))>1e-5 ...
+    || max(abs(oo_ramsey_policy_initval_AR2.dr.ys-oo_.dr.ys))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_initval_AR2.dr.ghx-oo_.dr.ghx)))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_initval_AR2.dr.ghu-oo_.dr.ghu)))>1e-5 ...
+    || max(max(abs(oo_ramsey_policy_initval_AR2.dr.Gy-oo_.dr.Gy)))>1e-5 ...
+    || max(abs((oo_ramsey_policy_initval_AR2.planner_objective_value-oo_.planner_objective_value)))>1e-5
+    error('Running stoch_simul after ramsey_policy leads to inconsistent results')
+end
\ No newline at end of file
diff --git a/tests/optimal_policy/Ramsey/ramsey_ex_wrong_ss_file.mod b/tests/optimal_policy/Ramsey/ramsey_ex_wrong_ss_file.mod
new file mode 100644
index 000000000..3b44bad17
--- /dev/null
+++ b/tests/optimal_policy/Ramsey/ramsey_ex_wrong_ss_file.mod
@@ -0,0 +1,41 @@
+/* Mod file tests whether the Ramsey command correctly filters out an incorrect steady state 
+ * The example is taken from Juillard, Michel (2011): User manual for optimal policy package, 
+ * MONFISPOL FP7 project SSH-225149, Deliverable 1.1.2
+*/
+
+
+var pai, c, n, r, a;
+varexo u;
+parameters beta, rho, epsilon, omega, phi, gamma;
+
+beta=0.99;
+gamma=3;
+omega=17;
+epsilon=8;
+phi=1;
+rho=0.95;
+
+model;
+a = rho*a(-1)+u;
+1/c = beta*r/(c(+1)*pai(+1));
+pai*(pai-1)/c = beta*pai(+1)*(pai(+1)-1)/c(+1)+epsilon*phi*n^(gamma+1)/omega -exp(a)*n*(epsilon-1)/(omega*c);
+exp(a)*n = c+(omega/2)*(pai-1)^2;
+end;
+
+initval;
+r=1;
+end;
+
+steady_state_model;
+a = 1;
+pai = beta*r;
+c = find_c(0.96,pai,beta,epsilon,phi,gamma,omega);
+n = c+(omega/2)*(pai-1)^2;
+end;
+
+shocks;
+var u; stderr 0.008;
+end;
+
+planner_objective(ln(c)-phi*((n^(1+gamma))/(1+gamma)));
+ramsey_policy(planner_discount=0.99,order=1,instruments=(r));
\ No newline at end of file