diff --git a/matlab/k_order_pert.m b/matlab/k_order_pert.m index c6a9d38b0..57d25bf50 100644 --- a/matlab/k_order_pert.m +++ b/matlab/k_order_pert.m @@ -39,6 +39,20 @@ switch(order) return; end dr.g_1 = g_1; + + dr.ghx = dr.g_1(:,1:nspred); + dr.ghu = dr.g_1(:,nspred+1:end); + + if options.loglinear == 1 + k = find(dr.kstate(:,2) <= M.maximum_endo_lag+1); + klag = dr.kstate(k,[1 2]); + k1 = dr.order_var; + + dr.ghx = repmat(1./dr.ys(k1),1,size(dr.ghx,2)).*dr.ghx.* ... + repmat(dr.ys(k1(klag(:,1)))',size(dr.ghx,1),1); + dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu; + end + case 2 [err, g_0, g_1, g_2] = k_order_perturbation(dr,M,options); if err @@ -48,8 +62,43 @@ switch(order) dr.g_0 = g_0; dr.g_1 = g_1; dr.g_2 = g_2; + + dr.ghx = dr.g_1(:,1:nspred); + dr.ghu = dr.g_1(:,nspred+1:end); + + dr.ghs2 = 2*g_0; + s0 = 0; + s1 = 0; + ghxx=zeros(endo_nbr, nspred^2); + ghxu=zeros(endo_nbr, nspred*exo_nbr); + ghuu=zeros(endo_nbr, exo_nbr^2); + for i=1:size(g_2,2) + if s0 < nspred && s1 < nspred + ghxx(:,s0*nspred+s1+1) = 2*g_2(:,i); + if s1 > s0 + ghxx(:,s1*nspred+s0+1) = 2*g_2(:,i); + end + elseif s0 < nspred && s1 < nspred+exo_nbr + ghxu(:,(s0*exo_nbr+s1-nspred+1)) = 2*g_2(:,i); + elseif s0 < nspred+exo_nbr && s1 < nspred+exo_nbr + ghuu(:,(s0-nspred)*exo_nbr+s1-nspred +1) = 2*g_2(:,i); + if s1 > s0 + ghuu(:,(s1-nspred)*exo_nbr+s0-nspred+1) = 2*g_2(:,i); + end + else + error('dr1:k_order_perturbation:g_2','Unaccounted columns in g_2'); + end + s1 = s1+1; + if s1 == nspred+exo_nbr + s0 = s0+1; + s1 = s0; + end + end % for loop + dr.ghxx = ghxx; + dr.ghxu = ghxu; + dr.ghuu = ghuu; case 3 - if options.pruning + if options.pruning [err, g_0, g_1, g_2, g_3, derivs] = k_order_perturbation(dr, ... M,options); if err @@ -84,58 +133,7 @@ switch(order) error('order > 3 isn''t implemented') end -if options.pruning - return -end -nspred = M.nspred; - -dr.ghx = dr.g_1(:,1:nspred); -dr.ghu = dr.g_1(:,nspred+1:end); - -if options.loglinear == 1 - k = find(dr.kstate(:,2) <= M.maximum_endo_lag+1); - klag = dr.kstate(k,[1 2]); - k1 = dr.order_var; - - dr.ghx = repmat(1./dr.ys(k1),1,size(dr.ghx,2)).*dr.ghx.* ... - repmat(dr.ys(k1(klag(:,1)))',size(dr.ghx,1),1); - dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu; -end - -if order > 1 - dr.ghs2 = 2*g_0; - s0 = 0; - s1 = 0; - ghxx=zeros(endo_nbr, nspred^2); - ghxu=zeros(endo_nbr, nspred*exo_nbr); - ghuu=zeros(endo_nbr, exo_nbr^2); - for i=1:size(g_2,2) - if s0 < nspred && s1 < nspred - ghxx(:,s0*nspred+s1+1) = 2*g_2(:,i); - if s1 > s0 - ghxx(:,s1*nspred+s0+1) = 2*g_2(:,i); - end - elseif s0 < nspred && s1 < nspred+exo_nbr - ghxu(:,(s0*exo_nbr+s1-nspred+1)) = 2*g_2(:,i); - elseif s0 < nspred+exo_nbr && s1 < nspred+exo_nbr - ghuu(:,(s0-nspred)*exo_nbr+s1-nspred +1) = 2*g_2(:,i); - if s1 > s0 - ghuu(:,(s1-nspred)*exo_nbr+s0-nspred+1) = 2*g_2(:,i); - end - else - error('dr1:k_order_perturbation:g_2','Unaccounted columns in g_2'); - end - s1 = s1+1; - if s1 == nspred+exo_nbr - s0 = s0+1; - s1 = s0; - end - end % for loop - dr.ghxx = ghxx; - dr.ghxu = ghxu; - dr.ghuu = ghuu; -end function y = unfold2(x,n) y=zeros(size(x,1),n*n);