diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m index 6b3ced3ad..3cf1a08c4 100644 --- a/matlab/ep/extended_path.m +++ b/matlab/ep/extended_path.m @@ -169,6 +169,9 @@ while (t3 @@ -152,7 +159,8 @@ if weight<1 solve_stochastic_perfect_foresight_model(endo_simul0,exo_simul0,pfm,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order); case 1 [flag,tmp] = ... - solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order); + solve_stochastic_perfect_foresight_model_1(endo_simul0,exxo_simul,options_.ep,pfm,options_.ep.stochastic.order,weight); + % solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order); end end end diff --git a/matlab/ep/setup_integration_nodes.m b/matlab/ep/setup_integration_nodes.m index da8e129b9..8907daac9 100644 --- a/matlab/ep/setup_integration_nodes.m +++ b/matlab/ep/setup_integration_nodes.m @@ -1,39 +1,39 @@ function [nodes,weights,nnodes] = setup_integration_nodes(EpOptions,pfm) if EpOptions.stochastic.order - % Compute weights and nodes for the stochastic version of the extended path. - switch EpOptions.IntegrationAlgorithm - case 'Tensor-Gaussian-Quadrature' - % Get the nodes and weights from a univariate Gauss-Hermite quadrature. - [nodes0,weights0] = gauss_hermite_weights_and_nodes(EpOptions.stochastic.quadrature.nodes); - % Replicate the univariate nodes for each innovation and dates, and, if needed, correlate them. - nodes0 = repmat(nodes0,1,pfm.number_of_shocks*pfm.stochastic_order)*kron(eye(pfm.stochastic_order),pfm.Omega); - % Put the nodes and weights in cells - for i=1:pfm.number_of_shocks - rr(i) = {nodes0(:,i)}; - ww(i) = {weights0}; + % Compute weights and nodes for the stochastic version of the extended path. + switch EpOptions.IntegrationAlgorithm + case 'Tensor-Gaussian-Quadrature' + % Get the nodes and weights from a univariate Gauss-Hermite quadrature. + [nodes0,weights0] = gauss_hermite_weights_and_nodes(EpOptions.stochastic.quadrature.nodes); + % Replicate the univariate nodes for each innovation and dates, and, if needed, correlate them. + nodes0 = repmat(nodes0,1,pfm.number_of_shocks*pfm.stochastic_order)*kron(eye(pfm.stochastic_order),pfm.Omega); + % Put the nodes and weights in cells + for i=1:pfm.number_of_shocks + rr(i) = {nodes0(:,i)}; + ww(i) = {weights0}; + end + % Build the tensorial grid + nodes = cartesian_product_of_sets(rr{:}); + weights = prod(cartesian_product_of_sets(ww{:}),2); + nnodes = length(weights); + case 'Stroud-Cubature-3' + [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,3,'Stroud') + nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; + weights = weights; + nnodes = length(weights); + case 'Stroud-Cubature-5' + [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,5,'Stroud') + nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; + weights = weights; + nnodes = length(weights); + case 'UT_2p+1' + p = pfm.number_of_shocks; + k = EpOptions.ut.k; + C = sqrt(pfm.number_of_shocks + k)*pfm.Omega'; + nodes = [zeros(1,p); -C; C]; + weights = [k/(p+k); (1/(2*(p+k)))*ones(2*p,1)]; + nnodes = 2*p+1; + otherwise + error('Stochastic extended path:: Unknown integration algorithm!') end - % Build the tensorial grid - nodes = cartesian_product_of_sets(rr{:}); - weights = prod(cartesian_product_of_sets(ww{:}),2); - nnodes = length(pfm.weights); - case 'Stroud-Cubature-3' - [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,3,'Stroud') - nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; - weights = weights; - nnodes = length(pfm.weights); - case 'Stroud-Cubature-5' - [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,5,'Stroud') - nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; - weights = weights; - nnodes = length(weights); - case 'UT_2p+1' - p = pfm.number_of_shocks; - k = EpOptions.ut.k; - C = sqrt(pfm.number_of_shocks + k)*pfm.Omega'; - nodes = [zeros(1,p); -C; C]; - weights = [k/(p+k); (1/(2*(p+k)))*ones(2*p,1)]; - nnodes = 2*p+1; - otherwise - error('Stochastic extended path:: Unknown integration algorithm!') end -end diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m index 7c7c3e5aa..b844f8ce8 100644 --- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m +++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m @@ -1,4 +1,4 @@ -function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,EpOptions,pfm,order) +function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,EpOptions,pfm,order,varargin) % Copyright (C) 2012-2013 Dynare Team % @@ -17,6 +17,11 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . + if nargin < 6 + homotopy_parameter = 1; + else + homotopy_parameter = varargin{1}; + end flag = 0; err = 0; stop = 0; @@ -159,7 +164,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo Y(i_cols_s,1); Y(i_cols_f,k1)]; end - [d1,jacobian] = dynamic_model(y,innovation, ... + [d1,jacobian] = dynamic_model(y,homotopy_parameter*innovation, ... params,steady_state,i+1); if i == 1 % in first period we don't keep track of @@ -182,7 +187,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo y = [Y(i_cols_p,1); Y(i_cols_s,j); Y(i_cols_f,j)]; - [d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1); + [d1,jacobian] = dynamic_model(y,homotopy_parameter*innovation,params,steady_state,i+1); i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; A1(i_rows,i_cols_A) = jacobian(:,i_cols_j); res(:,i,j) = d1; @@ -194,7 +199,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo y = [Y(i_cols_p,j); Y(i_cols_s,j); Y(i_cols_f,j)]; - [d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1); + [d1,jacobian] = dynamic_model(y,homotopy_parameter*innovation,params,steady_state,i+1); i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; A1(i_rows,i_cols_A) = jacobian(:,i_cols_j); res(:,i,j) = d1; @@ -242,17 +247,16 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo end if err < tolerance stop = 1; + flag = 0;% Convergency obtained. + endo_simul = reshape(Y(:,1),ny,periods+2);%Y(ny+(1:ny),1); if verbose + save ep_test_s1 exo_simul endo_simul Y fprintf('\n') ; disp([' Total time of simulation :' num2str(etime(clock,h1))]) ; fprintf('\n') ; disp([' Convergency obtained.']) ; fprintf('\n') ; end - flag = 0;% Convergency obtained. - endo_simul = reshape(Y(:,1),ny,periods+2);%Y(ny+(1:ny),1); - % figure;plot(Y(16:ny:(periods+2)*ny,:)) - % pause break end A2 = [nzA{:}]'; @@ -260,6 +264,8 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo if verbose disp(['solve_stochastic_foresight_model_1 encountered ' ... 'NaN']) + save ep_test_s2 exo_simul endo_simul + pause end flag = 1; return @@ -280,6 +286,8 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo disp(['WARNING : maximum number of iterations is reached (modify options_.simul.maxit).']) ; fprintf('\n') ; disp(sprintf('err: %f',err)); + save ep_test_s2 exo_simul endo_simul + pause end flag = 1;% more iterations are needed. endo_simul = 1;