diff --git a/matlab/discretionary_policy.m b/matlab/discretionary_policy.m
new file mode 100644
index 0000000000000000000000000000000000000000..eda6aa11dc4d80a881bd508de0ee10a5679bbe0f
--- /dev/null
+++ b/matlab/discretionary_policy.m
@@ -0,0 +1,40 @@
+function info = discretionary_policy(var_list)
+
+% Copyright (C) 2007-2011 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% 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_ oo_ M_
+
+options_.discretionary_policy = 1;
+oldoptions = options_;
+options_.order = 1;
+info = stoch_simul(var_list);
+
+if options_.noprint == 0
+    disp_steady_state(M_,oo_)
+    for i=M_.orig_endo_nbr:M_.endo_nbr
+        if strmatch('mult_',M_.endo_names(i,:))
+            disp(sprintf('%s \t\t %g',M_.endo_names(i,:), ...
+                         oo_.dr.ys(i)));
+        end
+    end
+end
+
+
+%oo_ = evaluate_planner_objective(oo_.dr,M_,oo_,options_);
+
+options_ = oldoptions;
\ No newline at end of file
diff --git a/matlab/discretionary_policy_1.m b/matlab/discretionary_policy_1.m
new file mode 100644
index 0000000000000000000000000000000000000000..fb9d226ec6da1173f08161c1e7308dcabbf4faa9
--- /dev/null
+++ b/matlab/discretionary_policy_1.m
@@ -0,0 +1,159 @@
+function [dr,ys,info]=discretionary_policy_1(oo_,Instruments)
+
+% Copyright (C) 2007-2011 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+global M_ options_
+persistent Hold
+
+options_ = set_default_option(options_,'noprint',0);
+options_ = set_default_option(options_,'qz_criterium',1.000001);
+options_ = set_default_option(options_,'solve_maxit',3000);
+options_ = set_default_option(options_,'discretion_tol',1e-7);
+
+% safeguard against issues like running ramsey policy first and then running discretion
+if isfield(M_,'orig_model')
+    orig_model = M_.orig_model;
+    M_.endo_nbr = orig_model.endo_nbr;
+    M_.endo_names = orig_model.endo_names;
+    M_.lead_lag_incidence = orig_model.lead_lag_incidence;
+    M_.maximum_lead = orig_model.maximum_lead;
+    M_.maximum_endo_lead = orig_model.maximum_endo_lead;
+    M_.maximum_lag = orig_model.maximum_lag;
+    M_.maximum_endo_lag = orig_model.maximum_endo_lag;
+else
+	M_.orig_model = M_;
+end
+
+beta = options_.planner_discount;
+
+exo_nbr = M_.exo_nbr;
+if isfield(M_,'orig_model')
+    orig_model = M_.orig_model;
+    endo_nbr = orig_model.endo_nbr;
+    endo_names = orig_model.endo_names;
+    lead_lag_incidence = orig_model.lead_lag_incidence;
+    MaxLead = orig_model.maximum_lead;
+    MaxLag = orig_model.maximum_lag;
+else
+	endo_names = M_.endo_names;
+	endo_nbr = M_.endo_nbr;
+	MaxLag=M_.maximum_lag;
+	MaxLead=M_.maximum_lead;
+    lead_lag_incidence = M_.lead_lag_incidence;
+end
+
+[U,Uy,W] = feval([M_.fname,'_objective_static'],zeros(endo_nbr,1),[], M_.params);
+if any(any(Uy~=0))
+    error(['discretionary_policy: the objective function must have zero ' ...
+           'first order derivatives'])
+end
+
+W=reshape(W,endo_nbr,endo_nbr);
+
+klen = MaxLag + MaxLead + 1;
+iyv=lead_lag_incidence';
+% Find the jacobian
+z = repmat(zeros(endo_nbr,1),1,klen);
+z = z(nonzeros(iyv)) ;
+it_ = MaxLag + 1 ;
+
+if exo_nbr == 0
+    oo_.exo_steady_state = [] ;
+end
+[junk,jacobia_] = feval([M_.fname '_dynamic'],z, [oo_.exo_simul ...
+                    oo_.exo_det_simul], M_.params, zeros(endo_nbr,1), it_);
+if any(junk~=0)
+    error(['discretionary_policy: the model must be written in deviation ' ...
+           'form and not have constant terms'])
+end
+
+eq_nbr= size(jacobia_,1);
+instr_nbr=endo_nbr-eq_nbr;
+
+instr_id=nan(instr_nbr,1);
+for j=1:instr_nbr
+	vj=deblank(Instruments(j,:));
+	vj_id=strmatch(vj,endo_names,'exact');
+	if ~isempty(vj_id)
+		instr_id(j)=vj_id;
+	else
+		error([mfilename,':: instrument ',vj,' not found'])
+	end
+end
+
+Indices={'lag','0','lead'};
+iter=1;
+for j=1:numel(Indices)
+    eval(['A',Indices{j},'=zeros(eq_nbr,endo_nbr);'])
+	if strcmp(Indices{j},'0')||(strcmp(Indices{j},'lag') && MaxLag)||(strcmp(Indices{j},'lead') && MaxLead)
+	    [junk,row,col]=find(lead_lag_incidence(iter,:));
+	    eval(['A',Indices{j},'(:,row)=jacobia_(:,col);'])
+		iter=iter+1;
+	end
+end
+B=jacobia_(:,nnz(iyv)+1:end);
+
+%%% MAIN ENGINE %%%
+qz_criterium = options_.qz_criterium;
+solve_maxit = options_.solve_maxit;
+discretion_tol = options_.discretion_tol;
+
+if ~isempty(Hold)
+	[H,G,info]=discretionary_policy_engine(Alag,A0,Alead,B,W,instr_id,beta,solve_maxit,discretion_tol,qz_criterium,Hold);
+else
+	[H,G,info]=discretionary_policy_engine(Alag,A0,Alead,B,W,instr_id,beta,solve_maxit,discretion_tol,qz_criterium);
+end
+
+if info
+    dr=[];
+    return
+else
+	Hold=H;
+	% Hold=[]; use this line if persistent command is not used.
+end
+% update the following elements
+
+LLI=lead_lag_incidence;
+LLI(MaxLag,:)=any(H);
+
+LLI=LLI';
+tmp=find(LLI);
+LLI(tmp)=1:numel(tmp);
+
+M_.lead_lag_incidence = LLI';
+
+% set the state
+dr=oo_.dr;
+dr.ys =zeros(endo_nbr,1);
+dr=set_state_space(dr,M_);
+order_var=dr.order_var;
+
+T=H(order_var,order_var);
+dr.ghu=G(order_var,:);
+Selection=any(T);
+dr.ghx=T(:,Selection);
+
+ys=NondistortionarySteadyState(M_);
+dr.ys=ys; % <--- dr.ys =zeros(NewEndo_nbr,1);
+
+function ys=NondistortionarySteadyState(M_)
+if exist([M_.fname,'_steadystate.m'],'file')
+	eval(['ys=',M_.fname,'_steadystate.m;'])
+else
+	ys=zeros(M_.endo_nbr,1);
+end
diff --git a/matlab/discretionary_policy_engine.m b/matlab/discretionary_policy_engine.m
new file mode 100644
index 0000000000000000000000000000000000000000..d8e109096f51505c2106e317beb397ff383df216
--- /dev/null
+++ b/matlab/discretionary_policy_engine.m
@@ -0,0 +1,297 @@
+function [H,G,retcode]=discretionary_policy_engine(AAlag,AA0,AAlead,BB,bigw,instr_id,beta,solve_maxit,discretion_tol,qz_criterium,H00,verbose)
+
+% Solves the discretionary problem for a model of the form:
+% AAlag*yy_{t-1}+AA0*yy_t+AAlead*yy_{t+1}+BB*e=0, with W the weight on the
+% variables in vector y_t and instr_id is the location of the instruments
+% in the yy_t vector.
+% We use the Dennis algorithm and so we need to re-write the model in the
+% form A0*y_t=A1*y_{t-1}+A2*y_{t+1}+A3*x_t+A4*x_{t+1}+A5*e_t, with W the
+% weight on the y_t vector and Q the weight on the x_t vector of
+% instruments.
+
+% Copyright (C) 2007-2011 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+if nargin<12
+    verbose=0;
+    if nargin<11
+        H00=[];
+        if nargin<10
+            qz_criterium=1.000001;
+            if nargin<9
+                discretion_tol=sqrt(eps);
+                if nargin<8
+                    solve_maxit=3000;
+                    if nargin<7
+                        beta=.99;
+                        if nargin<6
+                            error([mfilename,':: Insufficient number of input arguments'])
+                        elseif nargin>12
+                            error([mfilename,':: Number of input arguments cannot exceed 12'])
+                        end
+                    end
+                end
+            end
+        end
+    end
+end
+
+[A0,A1,A2,A3,A4,A5,W,Q,endo_nbr,exo_nbr,aux,endo_augm_id]=GetDennisMatrices(AAlag,AA0,AAlead,BB,bigw,instr_id);
+% aux is a logical index of the instruments which appear with lags in the
+% model. Their location in the state vector is instr_id(aux)
+% endo_augm_id is index (not logical) of locations of the augmented vector
+% of non-instrumental variables
+
+AuxiliaryVariables_nbr=sum(aux);
+H0=zeros(endo_nbr+AuxiliaryVariables_nbr);
+if ~isempty(H00)
+    H0(1:endo_nbr,1:endo_nbr)=H00;clear H00
+end
+
+H10=H0(endo_augm_id,endo_augm_id);
+F10=H0(instr_id,endo_augm_id);
+
+yy_nbr=numel(endo_augm_id);
+
+[H1,F1,retcode]=Iterate(H10,F10);
+
+switch retcode
+    case 2 % nan
+        retcode=3;
+        retcode(2)=10000;
+        if verbose
+            disp([mfilename,':: NAN elements in the solution'])
+        end
+    case 1% maxiter
+        if verbose
+            disp([mfilename,':: Maximum Number of Iterations reached'])
+        end
+    case 0
+        BadEig=max(abs(eig(H1)))>qz_criterium;
+        if BadEig
+            retcode=3;
+            retcode(2)=100*max(abs(eig(H1)));
+            if verbose
+                disp([mfilename,':: Some eigenvalues greater than qz_criterium, Model potentially unstable'])
+            end
+        end
+end
+
+if retcode(1)
+    H=[];
+    G=[];
+else
+    F2=-(Q+A3DPD*A3)\(A3DPD*A5);
+    H2=Dinv*(A5+A3*F2);
+    H=zeros(endo_nbr+AuxiliaryVariables_nbr);
+    G=zeros(endo_nbr+AuxiliaryVariables_nbr,exo_nbr);
+    H(endo_augm_id,endo_augm_id)=H1;
+    H(instr_id,endo_augm_id)=F1;
+    G(endo_augm_id,:)=H2;
+    G(instr_id,:)=F2;
+    
+    % Account for auxilliary variables
+    H(:,instr_id(aux))=H(:,end-(AuxiliaryVariables_nbr-1:-1:0));
+    H=H(1:endo_nbr,1:endo_nbr);
+    G=G(1:endo_nbr,:);
+end
+
+    function [H1,F1,retcode]=Iterate(H10,F10)
+        iter=0;
+        H1=H10;
+        F1=F10;
+        while 1
+            iter=iter+1;
+            P=SylvesterDoubling(W+beta*F1'*Q*F1,beta*H1',H1,discretion_tol,solve_maxit);
+            if any(any(isnan(P)))
+                P=SylvesterHessenbergSchur(W+beta*F1'*Q*F1,beta*H1',H1);
+                if any(any(isnan(P)))
+                    retcode=2;
+                    return
+                end
+            end
+            D=A0-A2*H1-A4*F1;
+            Dinv=D\eye(yy_nbr);
+            A3DPD=A3'*Dinv'*P*Dinv;
+            F1=-(Q+A3DPD*A3)\(A3DPD*A1);
+            H1=Dinv*(A1+A3*F1);
+            
+            [rcode,NQ]=CheckConvergence([H1;F1]-[H10;F10],iter,solve_maxit,discretion_tol);
+            if rcode
+                break
+            else
+                if verbose
+                    disp(NQ)
+                end
+            end
+            H10=H1;
+            F10=F1;
+        end
+        retcode=rcode-1;
+    end
+
+end
+
+function [rcode,NQ]=CheckConvergence(Q,iter,MaxIter,crit)
+
+NQ=max(max(abs(Q)));% norm(Q); seems too costly
+if isnan(NQ)
+    rcode=3;
+elseif iter>MaxIter;
+    rcode=2;
+elseif NQ<crit
+    rcode=1;
+else
+    rcode=0;
+end
+end
+
+function [A00,A11,A22,A33,A44,A55,WW,Q,endo_nbr,exo_nbr,aux,endo_augm_id]=GetDennisMatrices(AAlag,AA0,AAlead,BB,bigw,instr_id)
+[eq_nbr,endo_nbr]=size(AAlag);
+exo_nbr=size(BB,2);
+y=setdiff(1:endo_nbr,instr_id);
+instr_nbr=numel(instr_id);
+
+A0=AA0(:,y);
+A1=-AAlag(:,y);
+A2=-AAlead(:,y);
+A3=-AA0(:,instr_id);
+A4=-AAlead(:,instr_id);
+A5=-BB;
+W=bigw(y,y);
+Q=bigw(instr_id,instr_id);
+% Adjust for possible lags in instruments by creating auxiliary equations
+A6=-AAlag(:,instr_id);
+aux=any(A6);
+AuxiliaryVariables_nbr=sum(aux);
+ny=eq_nbr;
+m=eq_nbr+AuxiliaryVariables_nbr;
+A00=zeros(m);A00(1:ny,1:ny)=A0;A00(ny+1:end,ny+1:end)=eye(AuxiliaryVariables_nbr);
+A11=zeros(m);A11(1:ny,1:ny)=A1;A11(1:ny,ny+1:end)=A6(:,aux);
+A22=zeros(m);A22(1:ny,1:ny)=A2;
+A33=zeros(m,instr_nbr);A33(1:ny,1:end)=A3;A33(ny+1:end,aux)=eye(AuxiliaryVariables_nbr);
+A44=zeros(m,instr_nbr);A44(1:ny,1:end)=A4;
+A55=zeros(m,exo_nbr);A55(1:ny,1:end)=A5;
+WW=zeros(m);WW(1:ny,1:ny)=W;
+endo_augm_id=setdiff(1:endo_nbr+AuxiliaryVariables_nbr,instr_id);
+end
+
+function v= SylvesterDoubling (d,g,h,tol,maxit)
+
+% DOUBLES  Solves a Sylvester equation using doubling
+%
+%  [v,info] = doubles (g,d,h,tol,maxit)  uses a doubling algorithm
+%   to solve the  Sylvester equation v  = d + g v h
+
+v = d;
+for i =1:maxit,
+    vadd = g*v*h;
+    v = v+vadd;
+    if norm (vadd,1) <= (tol*norm(v,1))
+        break;
+    end
+    g = g*g;
+    h = h*h;
+end
+end
+
+function v = SylvesterHessenbergSchur(d,g,h)
+%
+% DSYLHS  Solves a discrete time sylvester equation	 using the
+% Hessenberg-Schur algorithm
+%
+% v = DSYLHS(g,d,h) computes the matrix v that satisfies the
+% discrete time sylvester equation
+%
+%           v = d + g'vh
+
+if size(g,1) >= size(h,1)
+    [u,gbarp] = hess(g');
+    [t,hbar] = schur(h);
+    [vbar] = sylvest_private(gbarp,u'*d*t,hbar,1e-15);
+    v = u*vbar*t';
+else
+    [u,gbar] = schur(g);
+    [t,hbarp] = hess(h');
+    [vbar] = sylvest_private(hbarp,t'*d'*u,gbar,1e-15);
+    v = u*vbar'*t';
+end
+end
+
+
+function v = sylvest_private(g,d,h,tol)
+%
+% SYLVEST  Solves a Sylvester equation
+%
+%  solves the Sylvester equation
+%         v = d + g v h
+%  for v where both g and h must be  upper block triangular.
+%  The output info is zero on a successful return.
+%  The input tol indicates when an element of g or h should be considered
+%  zero.
+
+[m,n] = size(d);
+v = zeros(m,n);
+w = eye(m);
+i = 1;
+temp = [];
+
+%First handle the i = 1 case outside the loop
+
+if i< n,
+    if abs(h(i+1,i)) < tol,
+        v(:,i)= (w - g*h(i,i))\d(:,i);
+        i = i+1;
+    else
+        A = [w-g*h(i,i)     (-g*h(i+1,i));...
+            -g*h(i,i+1)    w-g*h(i+1,i+1)];
+        C = [d(:,i); d(:,i+1)];
+        X = A\C;
+        v(:,i) = X(1:m,:);
+        v(:,i+1) = X(m+1:2*m,  :);
+        i = i+2;
+    end
+end
+
+%Handle the rest of the matrix with the possible exception of i=n
+
+while i<n,
+    b= i-1;
+    temp = [temp g*v(:,size(temp,2)+1:b)]; %#ok<AGROW>
+    if abs(h(i+1,i)) < tol,
+        v(:,i) = (w - g*h(i,i))\(d(:,i) + temp*h(1:b,i));
+        i = i+1;
+    else
+        A = [w - g*h(i,i)    (-g*h(i+1,i));   ...
+            -g*h(i,i+1)    w - g*h(i+1,i+1)];
+        C = [d(:,i) + temp*h(1:b,i);         ...
+            d(:,i+1) + temp*h(1:b,i+1)];
+        X = A\C;
+        v(:,i) = X(1:m,:);
+        v(:,i+1) = X(m+1:2*m, :);
+        i = i+2;
+    end
+end
+
+%Handle the i = n case if i=n was not in a 2-2 block
+
+if i==n,
+    b = i-1;
+    temp = [temp g*v(:,size(temp,2)+1:b)];
+    v(:,i) = (w-g*h(i,i))\(d(:,i) + temp*h(1:b,i));
+end
+end
\ No newline at end of file
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 013283ed701db5941529b82d9ce017e2267d8849..d82d6e3e1ef5fbef7e7590a83fddc42abaf87b8e 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -249,6 +249,9 @@ options_.initial_date.freq = 1;
 options_.initial_date.period = 1;
 options_.initial_date.subperiod = 0;
 
+% discretionary policy
+options_.discretionary_policy = 0;
+
 % SWZ SBVAR
 options_.ms.freq = 1;
 options_.ms.initial_subperiod = 1;
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 1fd8a07053ad74f8b893632af38bc55fc6488e96..7000c8a136693e8c83b6aab6ba934dd52f154a0d 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -62,6 +62,12 @@ check_model;
 
 if PI_PCL_solver
     [oo_.dr, info] = PCL_resol(oo_.steady_state,0);
+elseif options_.discretionary_policy
+    if ~options_.linear
+        error(['discretionary_policy solves only linear_quadratic ' ...
+               'problems']);
+    end
+    [oo_.dr,ys,info] = discretionary_policy_1(oo_,options_.instruments);
 else
     [oo_.dr, info] = resol(oo_.steady_state,0);
 end
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index e5fbb474ba3ae629d51ee21b680bd3b67bf4761b..8bde81382ecd2a75b1021156922044a0ad976297 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -216,6 +216,53 @@ RamseyPolicyStatement::writeOutput(ostream &output, const string &basename) cons
   output << "ramsey_policy(var_list_);\n";
 }
 
+DiscretionaryPolicyStatement::DiscretionaryPolicyStatement(const SymbolList &symbol_list_arg,
+							   const OptionsList &options_list_arg) :
+  symbol_list(symbol_list_arg),
+  options_list(options_list_arg)
+{
+}
+
+void
+DiscretionaryPolicyStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.discretionary_policy_present = true;
+
+  /* Fill in option_order of mod_file_struct
+     Since discretionary policy needs one further order of derivation (for example, for 1st order
+     approximation, it needs 2nd derivatives), we add 1 to the order declared by user */
+  OptionsList::num_options_t::const_iterator it = options_list.num_options.find("order");
+  if (it != options_list.num_options.end())
+    {
+      int order = atoi(it->second.c_str());
+      if (order > 1)
+        {
+          cerr << "ERROR: discretionary_policy: order > 1 is not yet implemented" << endl;
+          exit(EXIT_FAILURE);
+        }
+      mod_file_struct.order_option = max(mod_file_struct.order_option, order + 1);
+    }
+
+  // Fill in mod_file_struct.partial_information
+  it = options_list.num_options.find("partial_information");
+  if (it != options_list.num_options.end() && it->second == "1")
+    mod_file_struct.partial_information = true;
+
+  // Option k_order_solver (implicit when order >= 3)
+  it = options_list.num_options.find("k_order_solver");
+  if ((it != options_list.num_options.end() && it->second == "1")
+      || mod_file_struct.order_option >= 3)
+    mod_file_struct.k_order_solver = true;
+}
+
+void
+DiscretionaryPolicyStatement::writeOutput(ostream &output, const string &basename) const
+{
+  options_list.writeOutput(output);
+  symbol_list.writeOutput("var_list_", output);
+  output << "discretionary_policy(var_list_);\n";
+}
+
 EstimationStatement::EstimationStatement(const SymbolList &symbol_list_arg,
                                          const OptionsList &options_list_arg,
                                          const SymbolTable &symbol_table_arg) :
diff --git a/preprocessor/ComputingTasks.hh b/preprocessor/ComputingTasks.hh
index e265e95810496505784b666dadd2588fbac1c0b0..8be115492c4b757858e1fb1934f4a975a9b9e784 100644
--- a/preprocessor/ComputingTasks.hh
+++ b/preprocessor/ComputingTasks.hh
@@ -103,6 +103,18 @@ public:
   virtual void writeOutput(ostream &output, const string &basename) const;
 };
 
+class DiscretionaryPolicyStatement : public Statement
+{
+private:
+  const SymbolList symbol_list;
+  const OptionsList options_list;
+public:
+  DiscretionaryPolicyStatement(const SymbolList &symbol_list_arg,
+			       const OptionsList &options_list_arg);
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
 class RplotStatement : public Statement
 {
 private:
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 5a2a91e255a5eb0e9fd00bc4ed00fc3dc0771458..53b2a629c07089a6b152f54f899941b6625025f2 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -121,7 +121,7 @@ class ParsingDriver;
 %token RELATIVE_IRF REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE
 %token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO
 %token STDERR STEADY STOCH_SIMUL
-%token TEX RAMSEY_POLICY PLANNER_DISCOUNT
+%token TEX RAMSEY_POLICY PLANNER_DISCOUNT DISCRETIONARY_POLICY
 %token <string_val> TEX_NAME
 %token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL USEAUTOCORR
 %token VALUES VAR VAREXO VAREXO_DET VAROBS PREDETERMINED_VARIABLES
@@ -217,6 +217,7 @@ statement : parameters
           | model_info
           | planner_objective
           | ramsey_policy
+          | discretionary_policy
           | bvar_density
           | bvar_forecast
           | sbvar
@@ -1302,6 +1303,16 @@ ramsey_policy : RAMSEY_POLICY ';'
                 { driver.ramsey_policy(); }
               ;
 
+discretionary_policy : DISCRETIONARY_POLICY ';'
+                       { driver.discretionary_policy(); }
+                     | DISCRETIONARY_POLICY '(' ramsey_policy_options_list ')' ';'
+                       { driver.discretionary_policy(); }
+                     | DISCRETIONARY_POLICY symbol_list ';'
+                       { driver.discretionary_policy(); }
+                     | DISCRETIONARY_POLICY '(' ramsey_policy_options_list ')' symbol_list ';'
+                       { driver.discretionary_policy(); }
+                     ;
+
 ramsey_policy_options_list : ramsey_policy_options_list COMMA ramsey_policy_options
                            | ramsey_policy_options
                            ;
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 8e661a032e8aa50a154d2df11cd45cdcbd772079..1927bb50f77e16650ba2947ae298025d5f26b5be 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -130,6 +130,7 @@ string eofbuff;
 <INITIAL>Sigma_e {BEGIN DYNARE_STATEMENT; sigma_e = 1; return token::SIGMA_E;}
 <INITIAL>planner_objective {BEGIN DYNARE_STATEMENT; return token::PLANNER_OBJECTIVE;}
 <INITIAL>ramsey_policy {BEGIN DYNARE_STATEMENT; return token::RAMSEY_POLICY;}
+<INITIAL>discretionary_policy {BEGIN DYNARE_STATEMENT; return token::DISCRETIONARY_POLICY;}
 <INITIAL>identification {BEGIN DYNARE_STATEMENT; return token::IDENTIFICATION;}
 
 <INITIAL>bvar_density {BEGIN DYNARE_STATEMENT; return token::BVAR_DENSITY; }
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index 7fb56ce7d63878b629ab3a7ac3a3044f03e96ab5..ba9e818b2e797c1f2815d3cec255c64043f66377 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -1243,6 +1243,14 @@ ParsingDriver::ramsey_policy()
   options_list.clear();
 }
 
+void
+ParsingDriver::discretionary_policy()
+{
+  mod_file->addStatement(new DiscretionaryPolicyStatement(symbol_list, options_list));
+  symbol_list.clear();
+  options_list.clear();
+}
+
 void
 ParsingDriver::write_latex_dynamic_model()
 {
diff --git a/preprocessor/ParsingDriver.hh b/preprocessor/ParsingDriver.hh
index 50ea88f52e2ad5106a3f6e71095b0d3e0561565a..0491c44214416e88a1d03188553ff49471afab64 100644
--- a/preprocessor/ParsingDriver.hh
+++ b/preprocessor/ParsingDriver.hh
@@ -376,8 +376,10 @@ public:
   void begin_planner_objective();
   //! End a planner objective statement
   void end_planner_objective(expr_t expr);
-  //! ramsey policy statement
+  //! Ramsey policy statement
   void ramsey_policy();
+  //! Discretionary policy statement
+  void discretionary_policy();
   //! Adds a write_latex_dynamic_model statement
   void write_latex_dynamic_model();
   //! Adds a write_latex_static_model statement
diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc
index f8b4936bca47b4e89652247b4f8b921028e56018..ea34a236695beddc3a85dc3de7c1c0edc0a32a4d 100644
--- a/preprocessor/Statement.cc
+++ b/preprocessor/Statement.cc
@@ -29,6 +29,7 @@ ModFileStructure::ModFileStructure() :
   osr_params_present(false),
   optim_weights_present(false),
   ramsey_policy_present(false),
+  discretionary_policy_present(false),
   planner_objective_present(false),
   order_option(0),
   bvar_present(false),
diff --git a/preprocessor/Statement.hh b/preprocessor/Statement.hh
index a16f1646a0bada55b9a19eb3cb1a371d6c239d41..1d27815ebfff4b78a52fb0d4a0448910aad04c1b 100644
--- a/preprocessor/Statement.hh
+++ b/preprocessor/Statement.hh
@@ -50,6 +50,8 @@ public:
   bool optim_weights_present;
   //! Whether a ramsey_policy statement is present
   bool ramsey_policy_present;
+  //! Whether a discretionary_objective statement is present
+  bool discretionary_policy_present;
   //! Whether a planner_objective statement is present
   bool planner_objective_present;
   //! The value of the "order" option of stoch_simul, estimation, osr, ramsey_policy
diff --git a/tests/discretionary_policy/dennis_1.mod b/tests/discretionary_policy/dennis_1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..0643c07647bf3ff45fd571b091f11e21a5e9e6ea
--- /dev/null
+++ b/tests/discretionary_policy/dennis_1.mod
@@ -0,0 +1,37 @@
+var y i pi pi_c q;
+varexo g u e;
+
+parameters omega sigma beta kappa alpha;
+
+eta = 1;
+alpha = 0.4;
+beta = 0.99;
+theta = 0.75;
+sigma = 1;
+phi = 3;
+
+omega = 1+alpha*(2-alpha)*(sigma*eta-1);
+kappa = (1-theta)*(1-beta*theta)*(phi+sigma/omega)/theta;
+
+
+model(linear);
+y = y(+1) -(omega/sigma)*(i-pi(+1))+g;
+pi =  beta*pi(+1)+kappa*y+u;
+pi_c = pi+(alpha/(1-alpha))*(q-q(-1));
+q = q(+1)-(1-alpha)*(i-pi(+1))+(1-alpha)*e;
+i = 1.5*pi;
+end;
+
+shocks;
+var g; stderr 1;
+var u; stderr 1;
+var e; stderr 1;
+end;
+
+planner_objective pi_c^2 + y^2;
+
+discretionary_policy(instruments=(i),irf=0,qz_criterium=0.999999);
+
+return;
+
+ramsey_policy;
\ No newline at end of file