k_order_pert

 Compute decision rules using the k-order DLL from Dynare++

0001 function [dr,info] = k_order_pert(dr,M,options,oo)
0002 % Compute decision rules using the k-order DLL from Dynare++
0003 
0004 % Copyright (C) 2009-2012 Dynare Team
0005 %
0006 % This file is part of Dynare.
0007 %
0008 % Dynare is free software: you can redistribute it and/or modify
0009 % it under the terms of the GNU General Public License as published by
0010 % the Free Software Foundation, either version 3 of the License, or
0011 % (at your option) any later version.
0012 %
0013 % Dynare is distributed in the hope that it will be useful,
0014 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0015 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0016 % GNU General Public License for more details.
0017 %
0018 % You should have received a copy of the GNU General Public License
0019 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0020 
0021 info = 0;
0022 
0023 M.var_order_endo_names = M.endo_names(dr.order_var,:);
0024 
0025 order = options.order;
0026 
0027 switch(order)
0028   case 1
0029     [err, g_1] = k_order_perturbation(dr,M,options);
0030     mexErrCheck('k_order_perturbation', err);
0031     dr.g_1 = g_1;
0032   case 2
0033     [err, g_0, g_1, g_2] = k_order_perturbation(dr,M,options);
0034     mexErrCheck('k_order_perturbation', err);
0035     dr.g_0 = g_0;
0036     dr.g_1 = g_1;
0037     dr.g_2 = g_2;
0038   case 3
0039     [err, g_0, g_1, g_2, g_3] = k_order_perturbation(dr,M,options);
0040     mexErrCheck('k_order_perturbation', err);
0041     dr.g_0 = g_0;
0042     dr.g_1 = g_1;
0043     dr.g_2 = g_2;
0044     dr.g_3 = g_3;
0045   otherwise
0046     error('order > 3 isn''t implemented')
0047 end
0048 
0049 npred = dr.npred;
0050 
0051 dr.ghx = g_1(:,1:npred);
0052 dr.ghu = g_1(:,npred+1:end);
0053 
0054 if options.loglinear == 1
0055     k = find(dr.kstate(:,2) <= M.maximum_endo_lag+1);
0056     klag = dr.kstate(k,[1 2]);
0057     k1 = dr.order_var;
0058     
0059     dr.ghx = repmat(1./dr.ys(k1),1,size(dr.ghx,2)).*dr.ghx.* ...
0060              repmat(dr.ys(k1(klag(:,1)))',size(dr.ghx,1),1);
0061     dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu;
0062 end
0063 
0064 if order > 1
0065     dr.ghs2 = 2*g_0;
0066     endo_nbr = M.endo_nbr;
0067     exo_nbr = M.exo_nbr;
0068     s0 = 0;
0069     s1 = 0;
0070     ghxx=zeros(endo_nbr, npred^2);
0071     ghxu=zeros(endo_nbr, npred*exo_nbr);
0072     ghuu=zeros(endo_nbr, exo_nbr^2);
0073     for i=1:size(g_2,2)
0074         if s0 < npred && s1 < npred
0075             ghxx(:,s0*npred+s1+1) = 2*g_2(:,i);
0076             if s1 > s0
0077                 ghxx(:,s1*npred+s0+1) = 2*g_2(:,i);
0078             end
0079         elseif s0 < npred && s1 < npred+exo_nbr 
0080             ghxu(:,(s0*exo_nbr+s1-npred+1)) = 2*g_2(:,i);
0081         elseif s0 < npred+exo_nbr && s1 < npred+exo_nbr
0082             ghuu(:,(s0-npred)*exo_nbr+s1-npred +1) = 2*g_2(:,i);
0083             if s1 > s0
0084                 ghuu(:,(s1-npred)*exo_nbr+s0-npred+1) = 2*g_2(:,i);
0085             end
0086         else
0087             error('dr1:k_order_perturbation:g_2','Unaccounted columns in g_2');
0088         end
0089         s1 = s1+1;
0090         if s1 == npred+exo_nbr
0091             s0 = s0+1;
0092             s1 = s0; 
0093         end
0094     end % for loop
0095     dr.ghxx = ghxx;
0096     dr.ghxu = ghxu;
0097     dr.ghuu = ghuu;
0098 end
0099