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

Merge changes to extended path

time-shift
Michel Juillard 2014-05-12 09:43:49 +02:00
parent b5ca118bfb add594ab7e
commit c3efb214ef
9 changed files with 125 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
matlab/ep/extended_path.m
View File

@ -35,8 +35,11 @@ global M_ options_ oo_
options_.verbosity = options_.ep.verbosity;
verbosity = options_.ep.verbosity+options_.ep.debug;
% Set maximum number of iterations for the deterministic solver.
options_.simul.maxit = options_.ep.maxit;
% Prepare a structure needed by the matlab implementation of the perfect foresight model solver
pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_,'Tensor-Gaussian-Quadrature');
pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_);
exo_nbr = M_.exo_nbr;
periods = options_.periods;
@ -53,12 +56,10 @@ if isempty(initial_conditions)
end
end
% Set maximum number of iterations for the deterministic solver.
options_.simul.maxit = options_.ep.maxit;
% Set the number of periods for the perfect foresight model
periods = options_.ep.periods;
pfm.periods = options_.ep.periods;
pfm.periods = periods;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
% keep a copy of pfm.i_upd
@ -138,6 +139,17 @@ else
pfm.dr = [];
end
% number of nonzero derivatives
pfm.nnzA = M_.NNZDerivatives(1);
% setting up integration nodes if order > 0
if options_.ep.stochastic.order > 0
[nodes,weights,nnodes] = setup_integration_nodes(options_.ep,pfm);
pfm.nodes = nodes;
pfm.weights = weights;
pfm.nnodes = nnodes;
end
% Main loop.
while (t<sample_size)
if ~mod(t,10)
@ -168,6 +180,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
@ -185,7 +200,7 @@ while (t<sample_size)
solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
case 1
[flag,tmp] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
solve_stochastic_perfect_foresight_model_1(endo_simul_1,exo_simul_1,options_,pfm1,options_.ep.stochastic.order);
end
end
end
@ -263,7 +278,14 @@ while (t<sample_size)
if options_.ep.stochastic.order == 0
[flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
else
[flag,tmp] = solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.nodes,options_.ep.stochastic.order);
switch(options_.ep.stochastic.algo)
case 0
[flag,tmp] = ...
solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
case 1
[flag,tmp] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul_1,exo_simul_1,options_,pfm1,options_.ep.stochastic.order);
end
end
end
info_convergence = ~flag;

19
matlab/ep/homotopic_steps.m
View File

@ -33,12 +33,13 @@ exxo_simul = exo_simul0;
initial_step_length = step_length;
max_iter = 1000/step_length;
weight = initial_weight;
verbose = options_.ep.debug;
verbose = options_.ep.verbosity;
reduce_step_flag = 0;
if verbose
format long
disp('Entering homotopic_steps')
end
% (re)Set iter.
@ -59,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);
@ -69,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
@ -82,9 +87,10 @@ while weight<1
end
end
if info.convergence
%if d<stochastic_extended_path_depth
endo_simul0 = tmp;
%end
if options_.debug
save ep-test2 weight exo_simul0 endo_simul0
end
endo_simul0 = tmp;
jter = jter + 1;
if jter>3
if verbose
@ -153,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

39
matlab/ep/setup_integration_nodes.m Normal file
View File

@ -0,0 +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};
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
end

33
matlab/setup_stochastic_perfect_foresight_model_solver.m → matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
View File

@ -1,4 +1,4 @@
function pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,DynareOptions,DynareOutput,IntegrationMethod)
function pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,DynareOptions,DynareOutput)
% Copyright (C) 2013 Dynare Team
%
@ -69,34 +69,3 @@ pfm.verbose = DynareOptions.ep.verbosity;
pfm.maxit_ = DynareOptions.simul.maxit;
pfm.tolerance = DynareOptions.dynatol.f;
if nargin>3 && DynareOptions.ep.stochastic.order
% Compute weights and nodes for the stochastic version of the extended path.
switch IntegrationMethod
case 'Tensor-Gaussian-Quadrature'
% Get the nodes and weights from a univariate Gauss-Hermite quadrature.
[nodes,weights] = gauss_hermite_weights_and_nodes(DynareOptions.ep.stochastic.quadrature.nodes);
% Replicate the univariate nodes for each innovation and dates, and, if needed, correlate them.
nodes = repmat(nodes,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) = {nodes(:,i)};
ww(i) = {weights};
end
% Build the tensorial grid
pfm.nodes = cartesian_product_of_sets(rr{:});
pfm.weights = prod(cartesian_product_of_sets(ww{:}),2);
pfm.nnodes = length(pfm.weights);
case 'Stroud-Cubature-3'
[nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,3,'Stroud')
pfm.nodes = kron(eye(pfm.stochastic_order),transpose(Omega))*nodes;
pfm.weights = weights;
pfm.nnodes = length(pfm.weights);
case 'Stroud-Cubature-5'
[nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,5,'Stroud')
pfm.nodes = kron(eye(pfm.stochastic_order),transpose(Omega))*nodes;
pfm.weights = weights;
pfm.nnodes = length(weights);
otherwise
error('setup_stochastic_perfect_foresight_model_solver:: Unknown integration algorithm!')
end
end

70
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,pfm,nnodes,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;
@ -35,6 +40,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
hybrid_order = pfm.hybrid_order;
dr = pfm.dr;
[nodes,weights,nnodes] = setup_integration_nodes(EpOptions,pfm);
maxit = pfm.maxit_;
tolerance = pfm.tolerance;
@ -42,32 +48,16 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo
number_of_shocks = size(exo_simul,2);
[nodes,weights] = gauss_hermite_weights_and_nodes(nnodes);
% make sure that there is a node equal to zero
% and permute nodes and weights to have zero first
k = find(abs(nodes) < 1e-12);
k = find(sum(abs(nodes),2) < 1e-12);
if ~isempty(k)
nodes = [nodes(k); nodes(1:k-1); nodes(k+1:end)];
nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)];
weights = [weights(k); weights(1:k-1); weights(k+1:end)];
else
error('there is no nodes equal to zero')
end
if number_of_shocks>1
nodes = repmat(nodes,1,number_of_shocks)*chol(pfm.Sigma);
% to be fixed for Sigma ~= I
for i=1:number_of_shocks
rr(i) = {nodes(:,i)};
ww(i) = {weights};
end
nodes = cartesian_product_of_sets(rr{:});
weights = prod(cartesian_product_of_sets(ww{:}),2);
nnodes = nnodes^number_of_shocks;
else
nodes = nodes*sqrt(pfm.Sigma);
end
if hybrid_order > 0
if hybrid_order == 2
h_correction = 0.5*dr.ghs2(dr.inv_order_var);
@ -90,7 +80,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo
% The third row block is ny x nnodes^2
% and so on until size ny x nnodes^order
world_nbr = 1+(nnodes-1)*order;
Y = repmat(endo_simul(:),1,world_nbr);
Y = repmat(endo_simul(:),1,world_nbr);
% The columns of A map the elements of Y such that
% each block of Y with ny rows are unfolded column wise
@ -110,8 +100,8 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo
for i=2:periods
k = n1:n2;
for j=1:(1+(nnodes-1)*min(i-1,order))
i_upd_r(i1:i2) = (n1:n2)+(j-1)*ny*periods;
i_upd_y(i1:i2) = (n1:n2)+ny+(j-1)*ny*(periods+2);
i_upd_r(i1:i2) = k+(j-1)*ny*periods;
i_upd_y(i1:i2) = k+ny+(j-1)*ny*(periods+2);
i1 = i2+1;
i2 = i2+ny;
end
@ -174,7 +164,8 @@ 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,params,steady_state,i+1);
[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
% predetermined variables
@ -196,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;
@ -208,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;
@ -248,23 +239,41 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo
end
end
err = max(abs(res(i_upd_r)));
if verbose
[err1, k1] = max(abs(res));
[err2, k2] = max(abs(err1));
[err3, k3] = max(abs(err2));
disp([iter err k1(:,k2(:,:,k3),k3) k2(:,:,k3) k3])
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{:}]';
if any(isnan(A2(:,3))) || any(any(any(isnan(res))))
if verbose
disp(['solve_stochastic_foresight_model_1 encountered ' ...
'NaN'])
save ep_test_s2 exo_simul endo_simul
pause
end
flag = 1;
return
end
A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*world_nbr,dimension)];
if verbose
disp(sprintf('condest %g',condest(A)))
end
dy = -A\res(i_upd_r);
Y(i_upd_y) = Y(i_upd_y) + dy;
end
@ -276,6 +285,9 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo
fprintf('\n') ;
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;

1
preprocessor/DynamicModel.cc
View File

@ -4285,3 +4285,4 @@ DynamicModel::dynamicOnlyEquationsNbr() const
}

1
preprocessor/DynamicModel.hh
View File

@ -221,6 +221,7 @@ public:
void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order) const;
//! Writes file containing parameters derivatives
void writeParamsDerivativesFile(const string &basename) const;
//! Converts to static model (only the equations)
/*! It assumes that the static model given in argument has just been allocated */
void toStatic(StaticModel &static_model) const;

1
preprocessor/ModFile.cc
View File

@ -820,4 +820,3 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool no_log, b
cout << "done" << endl;
}

1
preprocessor/SteadyStateModel.cc
View File

@ -152,3 +152,4 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model
<< "end" << endl;
}