diff --git a/matlab/dyn_ramsey_static.m b/matlab/dyn_ramsey_static.m
index fd7b18d1233a59fb8c6d74086199d2754909f54b..564f2747b4f2630fd71d1ef52b59810177fa4d68 100644
--- a/matlab/dyn_ramsey_static.m
+++ b/matlab/dyn_ramsey_static.m
@@ -1,4 +1,4 @@
-function [steady_state,params,check] = dyn_ramsey_static(x,M,options_,oo)
+function [steady_state,params,check] = dyn_ramsey_static(ys_init,M,options_,oo)
 
 % function  [steady_state,params,check] = dyn_ramsey_static_(x)
 % Computes the static first order conditions for optimal policy
@@ -39,8 +39,8 @@ options_.steadystate.nocheck = 1; %disable checking because Lagrange multipliers
 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) && ~options_.steadystate_flag
-    steady_state = oo.steady_state;
+if check_static_model(ys_init,M,options_,oo) && ~options_.steadystate_flag
+    steady_state = ys_init;
     return
 elseif options_.steadystate_flag
     k_inst = [];
@@ -50,14 +50,14 @@ elseif options_.steadystate_flag
                                    M.endo_names,'exact')];
     end
     if inst_nbr == 1
-        inst_val = csolve(nl_func,oo.steady_state(k_inst),'',options_.solve_tolf,100);
+        inst_val = csolve(nl_func,ys_init(k_inst),'',options_.solve_tolf,100); %solve for instrument, using univariate solver, starting at initial value for instrument
     else
-        [inst_val,info1] = dynare_solve(nl_func,oo.steady_state(k_inst),0);
+        [inst_val,info1] = dynare_solve(nl_func,ys_init(k_inst),0); %solve for instrument, using multivariate solver, starting at initial value for instrument
     end
-    ys(k_inst) = inst_val;
+    ys_init(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,junk,steady_state] = nl_func(inst_val);
+    [xx,params,check] = evaluate_steady_state_file(ys_init,exo_ss,M,options_); %run steady state file again to update parameters
+    [junk,junk,steady_state] = nl_func(inst_val); %compute and return steady state
 else
     n_var = M.orig_endo_nbr;
     xx = oo.steady_state(1:n_var);
@@ -92,18 +92,19 @@ if options_.steadystate_flag
         k_inst = [k_inst; strmatch(instruments(i,:), ...
                                    endo_names,'exact')];
     end
-    oo.steady_state(k_inst) = x;
-    [x,params,check] = evaluate_steady_state_file(oo.steady_state,...
+    ys_init=zeros(size(oo.steady_state)); %create starting vector for steady state computation as only instrument value is handed over
+    ys_init(k_inst) = x; %set instrument, the only value required for steady state computation, to current value
+    [x,params,check] = evaluate_steady_state_file(ys_init,... returned x now has size endo_nbr as opposed to input size of n_instruments
                                                   [oo.exo_steady_state; ...
                                                   oo.exo_det_steady_state], ...
                                                   M,options_);
 end
 
-xx = zeros(endo_nbr,1);
-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)
+xx = zeros(endo_nbr,1); %initialize steady state vector
+xx(1:M.orig_endo_nbr) = x(1:M.orig_endo_nbr); %set values of original endogenous variables based on steady state file or initial value
+
+% setting steady state of auxiliary variables that depends on original endogenous variables
+if any([M.aux_vars.type] ~= 6) %auxiliary variables other than multipliers
     needs_set_auxiliary_variables = 1;
     fh = str2func([M.fname '_set_auxiliary_variables']);
     s_a_v_func = @(z) fh(z,... 
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 528f923797f91756626c2378d3c35443c4e1a316..3b295fb88982bbf8853fe84752fcea3254996208 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -29,6 +29,7 @@ MODFILES = \
 	optimal_policy/Ramsey/ramsey_ex_initval.mod \
 	optimal_policy/Ramsey/ramsey_ex.mod \
 	optimal_policy/Ramsey/ramsey_ex_initval_AR2.mod \
+	optimal_policy/Ramsey/ramsey_ex_aux.mod \
 	discretionary_policy/dennis_1.mod \
 	ramst_initval_file.mod \
 	ramst_normcdf_and_friends.mod \
diff --git a/tests/optimal_policy/Ramsey/ramsey_ex_aux.mod b/tests/optimal_policy/Ramsey/ramsey_ex_aux.mod
new file mode 100644
index 0000000000000000000000000000000000000000..bcf0669746306a7adac88af0564d87d1237993e4
--- /dev/null
+++ b/tests/optimal_policy/Ramsey/ramsey_ex_aux.mod
@@ -0,0 +1,52 @@
+/* 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, junk_endo_backward,junk_endo_forward,junk_expectation,junk_exo_backward,junk_exo_forward;
+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);
+junk_endo_backward=c(+2);
+junk_endo_forward=c(-2);
+junk_expectation=EXPECTATION(-1)(c(+1));
+junk_exo_backward=u(-2);
+junk_exo_forward=u(+2);
+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;
+junk_endo_backward=c;
+junk_endo_forward=c;
+junk_expectation=c;
+junk_exo_backward=0;
+junk_exo_forward=0;
+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));