stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

extended-path: changed homotopy so that integration nodes change 

too (only for for case 1)
time-shift
Michel Juillard 2014-03-16 18:49:07 +01:00
parent 8f06a3134d
commit f5ae90bcc8
4 changed files with 64 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
matlab/ep/extended_path.m
View File

@ -169,6 +169,9 @@ while (t<sample_size)
increase_periods = 0;
% Keep a copy of endo_simul_1
endo_simul = endo_simul_1;
if verbosity
save ep_test_1 endo_simul_1 exo_simul_1
end
while 1
if ~increase_periods
if bytecode_flag && ~options_.ep.stochastic.order

12
matlab/ep/homotopic_steps.m
View File

@ -60,6 +60,9 @@ while weight<1
else
flag = 1;
end
if options_.debug
save ep-test3 weight exo_simul0 endo_simul0
end
if flag
if ~options_.ep.stochastic.order
[flag,tmp,err] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
@ -70,7 +73,8 @@ while 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
@ -83,6 +87,9 @@ while weight<1
end
end
if info.convergence
if options_.debug
save ep-test2 weight exo_simul0 endo_simul0
end
endo_simul0 = tmp;
jter = jter + 1;
if jter>3
@ -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

70
matlab/ep/setup_integration_nodes.m
View File

@ -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

24
matlab/ep/solve_stochastic_perfect_foresight_model_1.m
View File

@ -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 <http://www.gnu.org/licenses/>.
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;