From cc5527c843a8bb8a2308f90a89602fc3a9757821 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Mon, 13 May 2013 19:23:38 +0200
Subject: [PATCH] Move filling of decision rules to respective orders

Previously, with pruning set to 1 and the option k_order_solver, the
fields like dr.ghu were not set in k_order_pert, even at order=1. But
stoch_simul tries to use them to display the decision rules. Similarly
for order 2. The respective setting of the decision rule fields is now
moved to the respective cases and always conducted at order 1 and
conducted at order 2. This avoids crashes if decision rules of
theoretical moments are requested.
---
 matlab/k_order_pert.m | 102 +++++++++++++++++++++---------------------
 1 file changed, 50 insertions(+), 52 deletions(-)

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);