(1) Added pruning algorithm for second order simulations (matlab code).

(2) Speed improvement in second order simulations (replaced call to matlab's kron function by call to
A_times_B_kronecker_C).

(3) Removed useless globals.

(4) Cosmetic changes and corrections on headers.

matlab/global_initialization.m

@@ -112,6 +112,7 @@ options_.minimal_solving_periods = 1;
 
 % Solution
 options_.order = 2;
+options_.pruning = 0;
 options_.solve_algo = 2;
 options_.linear = 0;
 options_.replic = 50;

matlab/simult.m

@@ -31,7 +31,6 @@ function y_=simult(ys, dr)
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 global M_ options_ oo_
-global  it_ means_
 
 order = options_.order;
 replic = options_.replic;
@@ -40,8 +39,6 @@ if replic == 0
 end
 seed = options_.simul_seed;
 
-it_ = M_.maximum_lag + 1 ;
-
 if replic > 1
     fname = [M_.fname,'_simul'];
     fh = fopen(fname,'w+');

matlab/simult_.m

@@ -1,23 +1,21 @@
 function y_=simult_(y0,dr,ex_,iorder)
-% function y_=simult_(y0,dr,ex_,iorder)
+% Simulates the model using a perturbation approach, given the path for the exogenous variables and the
-%
-% Simulates the model, given the path of exogenous variables and the
 % decision rules.
 %
 % INPUTS
-%    y0:       starting values
+%    y0       [double]   n*1 vector, initial value (n is the number of declared endogenous variables plus the number 
-%    dr:       structure of decisions rules for stochastic simulations
+%                        of auxilliary variables for lags and leads)
-%    ex_:      matrix of shocks
+%    dr       [struct]   matlab's structure where the reduced form solution of the model is stored.
-%    iorder=0: first-order approximation
+%    ex_      [double]   T*q matrix of innovations.
-%    iorder=1: second-order approximation
+%    iorder   [integer]  order of the taylor approximation.
 %
 % OUTPUTS
-%    y_:       stochastic simulations results
+%    y_       [double]   n*(T+1) time series for the endogenous variables.
 %
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2001-2007 Dynare Team
+% Copyright (C) 2001-2010 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -34,16 +32,15 @@ function y_=simult_(y0,dr,ex_,iorder)
 % 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_ it_
+global M_ options_
+
 iter = size(ex_,1);
-if ~isempty(dr.ghu)
+
-    nx = size(dr.ghu,2);
-end
 y_ = zeros(size(y0,1),iter+M_.maximum_lag);
 
-y_(:,1:M_.maximum_lag) = y0;
+y_(:,1) = y0;
 
-if options_.k_order_solver
+if options_.k_order_solver% Call dynare++ routines.
     options_.seed = 77;
     ex_ = [zeros(1,M_.exo_nbr); ex_];
     switch options_.order
@@ -64,46 +61,44 @@ if options_.k_order_solver
     end
     y_(dr.order_var,:) = y_;
 else
-    k1 = [M_.maximum_lag:-1:1];
     k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
     k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*M_.endo_nbr;
-    k3 = M_.lead_lag_incidence(1:M_.maximum_lag,:)';
+    switch iorder
-    k3 = find(k3(:));
+      case 1
-    k4 = dr.kstate(find(dr.kstate(:,2) < M_.maximum_lag+1),[1 2]);
+        if isempty(dr.ghu)
-    k4 = k4(:,1)+(M_.maximum_lag+1-k4(:,2))*M_.endo_nbr;
+            for i = 2:iter+M_.maximum_lag
-    if iorder == 1    
+                yhat = y_(dr.order_var(k2),i-1)-dr.ys(dr.order_var(k2));
-        if ~isempty(dr.ghu)
+                y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*yhat;
-            for i = M_.maximum_lag+1: iter+M_.maximum_lag
-                tempx1 = y_(dr.order_var,k1);
-                tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag);
-                tempx = tempx2(k2);
-                y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx+dr.ghu* ...
-                    ex_(i-M_.maximum_lag,:)';
-                k1 = k1+1;
             end
         else
-            for i = M_.maximum_lag+1: iter+M_.maximum_lag
+            for i = 2:iter+M_.maximum_lag
-                tempx1 = y_(dr.order_var,k1);
+                yhat = y_(dr.order_var(k2),i-1)-dr.ys(dr.order_var(k2));
-                tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag);
+                y_(dr.order_var,i) = dr.ys(dr.order_var) + dr.ghx*yhat + dr.ghu*ex_(i-1,:)';
-                tempx = tempx2(k2);
-                y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx;
-                k1 = k1+1;
             end
         end
-    elseif iorder == 2
+      case 2
-        for i = M_.maximum_lag+1: iter+M_.maximum_lag
+        constant = dr.ys(dr.order_var)+.5*dr.ghs2;
-            tempx1 = y_(dr.order_var,k1);
+        if options_.pruning
-            tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag);
+            y__ = y0;
-            tempx = tempx2(k2);
+            for i = 2:iter+M_.maximum_lag
-            tempu = ex_(i-M_.maximum_lag,:)';
+                yhat1 = y__(dr.order_var(k2))-dr.ys(dr.order_var(k2));
-            tempuu = kron(tempu,tempu);
+                yhat2 = y_(dr.order_var(k2),i-1)-dr.ys(dr.order_var(k2));
-            tempxx = kron(tempx,tempx);
+                epsilon = ex_(i-1,:)';
-            tempxu = kron(tempx,tempu);
+                y_(dr.order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
-            y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghs2/2+dr.ghx*tempx+ ...
+                    + A_times_B_kronecker_C(.5*dr.ghxx,yhat1) ...
-                dr.ghu*tempu+0.5*(dr.ghxx*tempxx+dr.ghuu*tempuu)+dr.ghxu*tempxu;
+                    + A_times_B_kronecker_C(.5*dr.ghuu,epsilon) ...
-            k1 = k1+1;
+                    + A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon);
+                y__(dr.order_var) = dr.ys(dr.order_var) + dr.ghx*yhat1 + dr.ghu*epsilon;
+            end
+        else
+            for i = 2:iter+M_.maximum_lag
+                yhat = y_(dr.order_var(k2),i-1)-dr.ys(dr.order_var(k2));
+                epsilon = ex_(i-1,:)';
+                y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
+                    + A_times_B_kronecker_C(.5*dr.ghxx,yhat) ...
+                    + A_times_B_kronecker_C(.5*dr.ghuu,epsilon) ...
+                    + A_times_B_kronecker_C(dr.ghxu,yhat,epsilon);
+            end
         end
     end
-end
+end
-
-% MJ 08/30/02 corrected bug at order 2