From 26f2b301b0d4777dc7127e3adc10a53e7994a6de Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Mon, 12 May 2014 10:18:43 +0200 Subject: [PATCH] make extended path algorithm 1 as a self contained problem usable by dynare_solve --- matlab/ep/ep_problem_2.m | 194 +++++++++ ...lve_stochastic_perfect_foresight_model_1.m | 398 ++++++------------ matlab/ep/stroud_judd_7.5.8.m | 9 + 3 files changed, 329 insertions(+), 272 deletions(-) create mode 100644 matlab/ep/ep_problem_2.m create mode 100644 matlab/ep/stroud_judd_7.5.8.m diff --git a/matlab/ep/ep_problem_2.m b/matlab/ep/ep_problem_2.m new file mode 100644 index 000000000..f646545a5 --- /dev/null +++ b/matlab/ep/ep_problem_2.m @@ -0,0 +1,194 @@ +function [res,A,info] = ep_problem_2(y,x,pm) + +info = 0; +res = []; +A = []; + +dynamic_model = pm.dynamic_model; +ny = pm.ny; +params = pm.params; +steady_state = pm.steady_state; +order = pm.order; +nodes = pm.nodes; +nnodes = pm.nnodes; +weights = pm.weights; +h_correction = pm.h_correction; +dimension = pm.dimension; +world_nbr = pm.world_nbr; +nnzA = pm.nnzA; +periods = pm.periods; +i_rows = pm.i_rows'; +i_cols = pm.i_cols; +nyp = pm.nyp; +nyf = pm.nyf; +hybrid_order = pm.hybrid_order; +i_cols_1 = pm.i_cols_1; +i_cols_j = pm.i_cols_j; +icA = pm.icA; +i_cols_T = pm.i_cols_T; + +i_cols_p = i_cols(1:nyp); +i_cols_s = i_cols(nyp+(1:ny)); +i_cols_f = i_cols(nyp+ny+(1:nyf)); +i_cols_A = i_cols; +i_cols_Ap0 = i_cols_p; +i_cols_As = i_cols_s; +i_cols_Af0 = i_cols_f - ny; +i_hc = i_cols_f - 2*ny; + +nzA = cell(periods,world_nbr); +res = zeros(ny,periods,world_nbr); +Y = zeros(ny*(periods+2),world_nbr); +Y(1:ny,1) = pm.y0; +Y(end-ny+1:end,:) = repmat(steady_state,1,world_nbr); +Y(pm.i_upd_y) = y; +offset_r0 = 0; +for i = 1:order+1 + i_w_p = 1; + for j = 1:(1+(nnodes-1)*(i-1)) + innovation = x; + if i <= order && j == 1 + % first world, integrating future shocks + if nargin > 1 + A1 = sparse([],[],[],i*ny,dimension,nnzA*world_nbr); + end + for k=1:nnodes + if nargin > 1 + if i == 2 + i_cols_Ap = i_cols_Ap0; + elseif i > 2 + i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes- ... + 1)*(i-2)*(i-3)/2); + end + if k == 1 + i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ... + 2); + else + i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ... + 2+(i-1)*(nnodes-1) ... + + k - 1); + end + end + if i > 1 + innovation(i+1,:) = nodes(k,:); + end + if k == 1 + k1 = 1; + else + k1 = (nnodes-1)*(i-1)+k; + end + if hybrid_order == 2 && (k > 1 || i == order) + z = [Y(i_cols_p,1); + Y(i_cols_s,1); + Y(i_cols_f,k1)+h_correction(i_hc)]; + else + z = [Y(i_cols_p,1); + Y(i_cols_s,1); + Y(i_cols_f,k1)]; + end + if nargin > 1 + [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1); + if i == 1 + % in first period we don't keep track of + % predetermined variables + i_cols_A = [i_cols_As - ny; i_cols_Af]; + A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1); + else + i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; + A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j); + end + else + d1 = dynamic_model(z,innovation,params,steady_state,i+1); + end + res(:,i,1) = res(:,i,1)+weights(k)*d1; + end + if nargin > 1 + [ir,ic,v] = find(A1); + nzA{i,j} = [ir,ic,v]'; + end + elseif j > 1 + (nnodes-1)*(i-2) + % new world, using previous state of world 1 and hit + % by shocks from nodes + if nargin > 1 + i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2); + i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); + end + k = j - (nnodes-1)*(i-2); + innovation(i+1,:) = nodes(k,:); + z = [Y(i_cols_p,1); + Y(i_cols_s,j); + Y(i_cols_f,j)]; + if nargin > 1 + [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1); + i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; + [ir,ic,v] = find(jacobian(:,i_cols_j)); + nzA{i,j} = [i_rows(ir),i_cols_A(ic), v]'; + else + d1 = dynamic_model(z,innovation,params,steady_state,i+1); + end + res(:,i,j) = d1; + if nargin > 1 + i_cols_Af = i_cols_Af + ny; + end + else + % existing worlds other than 1 + if nargin > 1 + i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2+j-1); + i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); + end + z = [Y(i_cols_p,j); + Y(i_cols_s,j); + Y(i_cols_f,j)]; + if nargin > 1 + [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1); + i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; + [ir,ic,v] = find(jacobian(:,i_cols_j)); + nzA{i,j} = [i_rows(ir),i_cols_A(ic),v]'; + else + d1 = dynamic_model(z,innovation,params,steady_state,i+1); + end + res(:,i,j) = d1; + i_cols_Af = i_cols_Af + ny; + end + i_rows = i_rows + ny; + if mod(j,nnodes) == 0 + i_w_p = i_w_p + 1; + end + if nargin > 1 && i > 1 + i_cols_As = i_cols_As + ny; + end + offset_r0 = offset_r0 + ny; + end + i_cols_p = i_cols_p + ny; + i_cols_s = i_cols_s + ny; + i_cols_f = i_cols_f + ny; +end +for j=1:world_nbr + i_rows_y = i_cols+(order+1)*ny; + offset_c = ny*(order+(nnodes-1)*(order-1)*order/2+j-1); + offset_r = offset_r0+(j-1)*ny; + for i=order+2:periods + if nargin > 1 + [d1,jacobian] = dynamic_model(Y(i_rows_y,j),x,params, ... + steady_state,i+1); + if i < periods + [ir,ic,v] = find(jacobian(:,i_cols_j)); + else + [ir,ic,v] = find(jacobian(:,i_cols_T)); + end + nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]'; + else + d1 = dynamic_model(Y(i_rows_y,j),x,params, ... + steady_state,i+1); + end + res(:,i,j) = d1; + i_rows_y = i_rows_y + ny; + offset_c = offset_c + world_nbr*ny; + offset_r = offset_r + world_nbr*ny; + end +end +if nargin > 1 + iA = [nzA{:}]'; + A = sparse(iA(:,1),iA(:,2),iA(:,3),dimension,dimension); +end +res = res(pm.i_upd_r); \ No newline at end of file diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m index b844f8ce8..ddebc0dbe 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,varargin) +function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,Options,pfm,order,varargin) % Copyright (C) 2012-2013 Dynare Team % @@ -17,281 +17,135 @@ 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}; +global options_ + +if nargin < 6 + homotopy_parameter = 1; +else + homotopy_parameter = varargin{1}; +end +flag = 0; +err = 0; +stop = 0; + +EpOptions = Options.ep; + +params = pfm.params; +steady_state = pfm.steady_state; +ny = pfm.ny; +periods = pfm.periods; +dynamic_model = pfm.dynamic_model; +lead_lag_incidence = pfm.lead_lag_incidence; +nyp = pfm.nyp; +nyf = pfm.nyf; +i_cols_1 = pfm.i_cols_1; +i_cols_A1 = pfm.i_cols_A1; +i_cols_j = pfm.i_cols_j; +i_cols_T = nonzeros(lead_lag_incidence(1:2,:)'); +hybrid_order = pfm.hybrid_order; +dr = pfm.dr; +nodes = pfm.nodes; +weights = pfm.weights; +nnodes = pfm.nnodes; + +maxit = pfm.maxit_; +tolerance = pfm.tolerance; +verbose = pfm.verbose; + +number_of_shocks = size(exo_simul,2); + +% make sure that there is a node equal to zero +% and permute nodes and weights to have zero first +k = find(sum(abs(nodes),2) < 1e-12); +if ~isempty(k) + 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 hybrid_order > 0 + if hybrid_order == 2 + h_correction = 0.5*dr.ghs2(dr.inv_order_var); end - flag = 0; - err = 0; - stop = 0; +else + h_correction = 0; +end - params = pfm.params; - steady_state = pfm.steady_state; - ny = pfm.ny; - periods = pfm.periods; - dynamic_model = pfm.dynamic_model; - lead_lag_incidence = pfm.lead_lag_incidence; - nyp = pfm.nyp; - nyf = pfm.nyf; - i_cols_1 = pfm.i_cols_1; - i_cols_A1 = pfm.i_cols_A1; - i_cols_j = pfm.i_cols_j; - 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; - verbose = pfm.verbose; +if verbose + disp ([' -----------------------------------------------------']); + disp (['MODEL SIMULATION :']); + fprintf('\n'); +end - number_of_shocks = size(exo_simul,2); +% Each column of Y represents a different world +% The upper right cells are unused +% The first row block is ny x 1 +% The second row block is ny x nnodes +% 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 = endo_simul(:,2:end-1); +Y = repmat(Y,1,world_nbr); +pfm.y0 = endo_simul(:,1); - % make sure that there is a node equal to zero - % and permute nodes and weights to have zero first - k = find(sum(abs(nodes),2) < 1e-12); - if ~isempty(k) - 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 hybrid_order > 0 - if hybrid_order == 2 - h_correction = 0.5*dr.ghs2(dr.inv_order_var); - end +% The columns of A map the elements of Y such that +% each block of Y with ny rows are unfolded column wise +% number of blocks +block_nbr = (order+(nnodes-1)*(order-1)*order/2+(periods-order)*world_nbr); +% dimension of the problem +dimension = ny*block_nbr; +pfm.block_nbr = block_nbr; +pfm.dimension = dimension; +if order == 0 + i_upd_r = (1:ny*periods); + i_upd_y = i_upd_r + ny; +else + i_upd_r = zeros(dimension,1); + i_upd_y = i_upd_r; + i_upd_r(1:ny) = (1:ny); + i_upd_y(1:ny) = ny+(1:ny); + i1 = ny+1; + i2 = 2*ny; + n1 = ny+1; + n2 = 2*ny; + for i=2:periods + k = n1:n2; + for j=1:(1+(nnodes-1)*min(i-1,order)) + 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 + n1 = n2+1; + n2 = n2+ny; end +end +icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ... + find(lead_lag_incidence(3,:))+2*world_nbr*ny]'; +h1 = clock; +pfm.order = order; +pfm.world_nbr = world_nbr; +pfm.nodes = nodes; +pfm.nnodes = nnodes; +pfm.weights = weights; +pfm.h_correction = h_correction; +pfm.i_rows = 1:ny; +i_cols = find(lead_lag_incidence'); +pfm.i_cols = i_cols; +pfm.nyp = nyp; +pfm.nyf = nyf; +pfm.hybrid_order = hybrid_order; +pfm.i_cols_1 = i_cols_1; +pfm.i_cols_h = i_cols_j; +pfm.icA = icA; +pfm.i_cols_T = i_cols_T; +pfm.i_upd_r = i_upd_r; +pfm.i_upd_y = i_upd_y; - if verbose - disp ([' -----------------------------------------------------']); - disp (['MODEL SIMULATION :']); - fprintf('\n'); - end - z = endo_simul(find(lead_lag_incidence')); - [d1,jacobian] = dynamic_model(z,exo_simul,params,steady_state,2); - - % Each column of Y represents a different world - % The upper right cells are unused - % The first row block is ny x 1 - % The second row block is ny x nnodes - % 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); - - % The columns of A map the elements of Y such that - % each block of Y with ny rows are unfolded column wise - dimension = ny*(order+(nnodes-1)*(order-1)*order/2+(periods-order)*world_nbr); - if order == 0 - i_upd_r = (1:ny*periods); - i_upd_y = i_upd_r + ny; - else - i_upd_r = zeros(dimension,1); - i_upd_y = i_upd_r; - i_upd_r(1:ny) = (1:ny); - i_upd_y(1:ny) = ny+(1:ny); - i1 = ny+1; - i2 = 2*ny; - n1 = ny+1; - n2 = 2*ny; - for i=2:periods - k = n1:n2; - for j=1:(1+(nnodes-1)*min(i-1,order)) - 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 - n1 = n2+1; - n2 = n2+ny; - end - end - icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ... - find(lead_lag_incidence(3,:))+2*world_nbr*ny]'; - h1 = clock; - for iter = 1:maxit - h2 = clock; - A1 = sparse([],[],[],ny*(order+(nnodes-1)*(order-1)*order/2),dimension,(order+1)*world_nbr*nnz(jacobian)); - res = zeros(ny,periods,world_nbr); - i_rows = 1:ny; - i_cols = find(lead_lag_incidence'); - i_cols_p = i_cols(1:nyp); - i_cols_s = i_cols(nyp+(1:ny)); - i_cols_f = i_cols(nyp+ny+(1:nyf)); - i_cols_A = i_cols; - i_cols_Ap0 = i_cols_p; - i_cols_As = i_cols_s; - i_cols_Af0 = i_cols_f - ny; - i_hc = i_cols_f - 2*ny; - for i = 1:order+1 - i_w_p = 1; - for j = 1:(1+(nnodes-1)*(i-1)) - innovation = exo_simul; - if i <= order && j == 1 - % first world, integrating future shocks - for k=1:nnodes - if i == 2 - i_cols_Ap = i_cols_Ap0; - elseif i > 2 - i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes- ... - 1)*(i-2)*(i-3)/2); - end - if k == 1 - i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ... - 2); - else - i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ... - 2+(i-1)*(nnodes-1) ... - + k - 1); - end - if i > 1 - innovation(i+1,:) = nodes(k,:); - end - if k == 1 - k1 = 1; - else - k1 = (nnodes-1)*(i-1)+k; - end - if hybrid_order == 2 && (k > 1 || i == order) - y = [Y(i_cols_p,1); - Y(i_cols_s,1); - Y(i_cols_f,k1)+h_correction(i_hc)]; - else - y = [Y(i_cols_p,1); - Y(i_cols_s,1); - Y(i_cols_f,k1)]; - end - [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 - i_cols_A = [i_cols_As - ny; i_cols_Af]; - A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1); - else - i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; - A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j); - end - res(:,i,1) = res(:,i,1)+weights(k)*d1; - end - elseif j > 1 + (nnodes-1)*(i-2) - % new world, using previous state of world 1 and hit - % by shocks from nodes - i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2); - i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); - k = j - (nnodes-1)*(i-2); - innovation(i+1,:) = nodes(k,:); - y = [Y(i_cols_p,1); - Y(i_cols_s,j); - Y(i_cols_f,j)]; - [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; - i_cols_Af = i_cols_Af + ny; - else - % existing worlds other than 1 - i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2+j-1); - i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); - y = [Y(i_cols_p,j); - Y(i_cols_s,j); - Y(i_cols_f,j)]; - [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; - i_cols_Af = i_cols_Af + ny; - end - i_rows = i_rows + ny; - if mod(j,nnodes) == 0 - i_w_p = i_w_p + 1; - end - if i > 1 - i_cols_As = i_cols_As + ny; - end - end - i_cols_p = i_cols_p + ny; - i_cols_s = i_cols_s + ny; - i_cols_f = i_cols_f + ny; - end - nzA = cell(periods,world_nbr); - parfor j=1:world_nbr - i_rows_y = find(lead_lag_incidence')+(order+1)*ny; - offset_c = ny*(order+(nnodes-1)*(order-1)*order/2+j-1); - offset_r = (j-1)*ny; - for i=order+2:periods - [d1,jacobian] = dynamic_model(Y(i_rows_y,j), ... - exo_simul,params, ... - steady_state,i+1); - if i == periods - [ir,ic,v] = find(jacobian(:,i_cols_T)); - else - [ir,ic,v] = find(jacobian(:,i_cols_j)); - end - nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]'; - res(:,i,j) = d1; - i_rows_y = i_rows_y + ny; - offset_c = offset_c + world_nbr*ny; - offset_r = offset_r + world_nbr*ny; - 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 - 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 - - if ~stop - if verbose - fprintf('\n') ; - disp([' Total time of simulation :' num2str(etime(clock,h1))]) ; - 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; - end - if verbose - disp (['-----------------------------------------------------']) ; - end \ No newline at end of file +options_.solve_algo = 9; +options_.steady.maxit = 100; +y = repmat(steady_state,block_nbr,1); +y = dynare_solve(@ep_problem_2,y,1,exo_simul,pfm); +endo_simul(:,2) = y(1:ny); \ No newline at end of file diff --git a/matlab/ep/stroud_judd_7.5.8.m b/matlab/ep/stroud_judd_7.5.8.m new file mode 100644 index 000000000..40c9bc08d --- /dev/null +++ b/matlab/ep/stroud_judd_7.5.8.m @@ -0,0 +1,9 @@ +function [X,w]=stroud_judd_7.5.8(d) + + E = eye(d); + X = cell(2*d,1); + m = 1; + for i=1:d + X{m} = E(:,i); + m = m+1; + X{m} = -E(:,i);