diff --git a/matlab/forcst2.m b/matlab/forcst2.m
index 41bf18f48..efa5d463b 100644
--- a/matlab/forcst2.m
+++ b/matlab/forcst2.m
@@ -1,5 +1,16 @@
 function yf=forcst2(y0,horizon,dr,n)
-
+% function yf=forcst2(y0,horizon,dr,n)
+% 
+% computes forecasts based on first order model solution, given shocks drawn from the shock distribution
+% Inputs:
+%   - y0        [endo_nbr by maximum_endo_lag]          matrix of starting values
+%   - dr        [structure]                         structure with Dynare decision rules
+%   - e         [horizon by exo_nbr]                matrix with shocks
+%   - n         [scalar]                            number of repetitions
+%
+% Outputs:
+%   - yf        [horizon+ykmin_ by endo_nbr by n]   array of forecasts
+% 
 % Copyright (C) 2008-2010 Dynare Team
 %
 % This file is part of Dynare.
@@ -24,14 +35,10 @@ endo_nbr = M_.endo_nbr;
 exo_nbr = M_.exo_nbr;
 ykmin_ = M_.maximum_endo_lag;
 
-order = options_.order;
-
 k1 = [ykmin_:-1:1];
 k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
 k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr;
 
-it_ = ykmin_ + 1 ;
-
 % eliminate shocks with 0 variance
 i_exo_var = setdiff([1:exo_nbr],find(diag(Sigma_e_) == 0));
 nxs = length(i_exo_var);
diff --git a/matlab/forcst2a.m b/matlab/forcst2a.m
index ecfb665af..16fcdc143 100644
--- a/matlab/forcst2a.m
+++ b/matlab/forcst2a.m
@@ -1,6 +1,15 @@
 function yf=forcst2a(y0,dr,e)
-
-% Copyright (C) 2008-2010 Dynare Team
+% function yf=forcst2a(y0,dr,e)
+% computes forecasts based on first order model solution, assuming the absence of shocks 
+% Inputs:
+%   - y0        [endo_nbr by maximum_endo_lag]          matrix of starting values
+%   - dr        [structure]                             structure with Dynare decision rules
+%   - e         [horizon by exo_nbr]                    matrix with shocks
+%
+% Outputs:
+%   - yf        [horizon+maximum_endo_lag,endo_nbr]               matrix of forecasts
+%
+% Copyright (C) 2008-2015 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -19,13 +28,10 @@ function yf=forcst2a(y0,dr,e)
 
 global M_ options_
 
-Sigma_e_ = M_.Sigma_e;
 endo_nbr = M_.endo_nbr;
-exo_nbr = M_.exo_nbr;
 ykmin_ = M_.maximum_endo_lag;
 
 horizon = size(e,1);
-order = options_.order;
 
 k1 = [ykmin_:-1:1];
 k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);