Added Stochastic Extended Path method (SEP).

matlab/ep/extended_path.m

@@ -79,7 +79,7 @@ make_ex_;
 make_y_;
 
 % Initialize the output array.
-time_series = NaN(M_.endo_nbr,sample_size+1);
+time_series = zeros(M_.endo_nbr,sample_size);
 
 % Set the covariance matrix of the structural innovations.
 variances = diag(M_.Sigma_e); 
@@ -106,6 +106,7 @@ if options_.ep.stochastic
     [r,w] = gauss_hermite_weights_and_nodes(options_.ep.number_of_nodes);
     switch options_.ep.stochastic
       case 1
+        if M_.exo_nbr>1
             rr = cell(1);
             ww = cell(1);
             for i=1:size(M_.Sigma_e,1)
@@ -114,10 +115,19 @@ if options_.ep.stochastic
             end
             rrr = cartesian_product_of_sets(rr{:});
             www = cartesian_product_of_sets(ww{:});
+        else
+            rrr = r;
+            www = w;
+        end
         www = prod(www,2);
+        nnn = length(www);    
       otherwise
         error(['Order ' int2str(options_.ep.stochastic) ' Stochastic Extended Path method is not implemented!'])
     end
+else
+    rrr = zeros(1,number_of_structural_innovations);
+    www = 1;
+    nnn = 1;
 end
 
 % Initializes some variables.
@@ -158,7 +168,9 @@ while (t<sample_size)
     shocks = oo_.ep.shocks(t,:);
     % Put it in oo_.exo_simul (second line).
     oo_.exo_simul(2,positive_var_indx) = shocks;
-    if options_.ep.init && t==1% Compute first order solution.
+    for s = 1:nnn
+        oo_.exo_simul(3,positive_var_indx) = rrr(s,:)*covariance_matrix_upper_cholesky;
+        if options_.ep.init && s==1% Compute first order solution. t==1 && 
             initial_path = simult_(initial_conditions,oo_.dr,oo_.exo_simul(2:end,:),1);
             if options_.ep.init==1
                 oo_.endo_simul(:,1:end-1) = initial_path(:,1:end-1);% Last column is the steady state.
@@ -278,20 +290,15 @@ while (t<sample_size)
         else
             oo_.endo_simul = tmp;
         end
-    if nargin==6
-        zlb_periods = find(oo_.endo_simul(zlb_idx,:)<=1+1e-12);
-        zlb_number_of_periods = length(zlb_periods);
-        if zlb_number_of_periods
-            count_zlb = [count_zlb ; [t, zlb_number_of_periods, zlb_periods(1) , zlb_periods(end)] ];
-        end
-    end 
         % Save results of the perfect foresight model solver.
-    time_series(:,t) = oo_.endo_simul(:,2);
-    save('simulated_paths.mat','time_series');
+        time_series(:,t) = time_series(:,t)+ www(s)*oo_.endo_simul(:,2);
+        %save('simulated_paths.mat','time_series');
         % Set initial condition for the nex round.
         %initial_conditions = oo_.endo_simul(:,2);
-    oo_.endo_simul = oo_.endo_simul(:,1:options_.periods+2);
+    end
+    %oo_.endo_simul = oo_.endo_simul(:,1:options_.periods+M_.maximum_endo_lag+M_.maximum_endo_lead);
     oo_.endo_simul(:,1:end-1) = oo_.endo_simul(:,2:end);
+    oo_.endo_simul(:,1) = time_series(:,t);
     oo_.endo_simul(:,end) = oo_.steady_state;
 end
 
@@ -302,3 +309,5 @@ else
         fprintf(back);
     end
 end
+
+oo_.endo_simul = oo_.steady_state;
\ No newline at end of file