From 48dad3e37aac222fff7416f1dc33f1d8c2424228 Mon Sep 17 00:00:00 2001 From: michel Date: Mon, 30 Nov 2009 19:54:45 +0000 Subject: [PATCH] 4.1: add preprocessor interface for k_order_solver rename use_k_order to k_order_solver in Matlab procedures fix tests git-svn-id: https://www.dynare.org/svn/dynare/trunk@3179 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/dr1.m | 2 +- matlab/global_initialization.m | 2 +- matlab/k_order_pert.m | 14 ++-- matlab/simult_.m | 90 ++++++++++++++---------- matlab/stoch_simul.m | 2 +- preprocessor/DynareBison.yy | 4 +- preprocessor/DynareFlex.ll | 1 + tests/k_order_perturbation/fs2000k2.mod | 2 +- tests/k_order_perturbation/fs2000k_1.mod | 32 +++++---- 9 files changed, 88 insertions(+), 61 deletions(-) diff --git a/matlab/dr1.m b/matlab/dr1.m index f413c0644..82c22690c 100644 --- a/matlab/dr1.m +++ b/matlab/dr1.m @@ -51,7 +51,7 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_) info = 0; - if options_.use_k_order; + if options_.k_order_solver; dr = set_state_space(dr,M_); [dr,info] = k_order_pert(dr,M_,options_,oo_); oo_.dr = dr; diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 7155eb0e1..d35ccb329 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -127,7 +127,7 @@ function global_initialization() options_.use_qzdiv = 0; end options_.aim_solver = 0; % i.e. by default do not use G.Anderson's AIM solver, use mjdgges instead - options_.use_k_order=0; % by default do not use k_order_perturbation but mjdgges + options_.k_order_solver=0; % by default do not use k_order_perturbation but mjdgges options_.partial_information = 0; options_.conditional_variance_decomposition = []; diff --git a/matlab/k_order_pert.m b/matlab/k_order_pert.m index b0f53c9ff..2094bec8a 100644 --- a/matlab/k_order_pert.m +++ b/matlab/k_order_pert.m @@ -10,12 +10,20 @@ function [dr,info] = k_order_pert(dr,M,options,oo) case 1 g_1 = k_order_perturbation(dr,0,M,options, oo , ['.' ... mexext]); + dr.g_1 = g_1; case 2 [g_0, g_1, g_2] = k_order_perturbation(dr,0,M,options, oo , ['.' ... mexext]); + dr.g_0 = g_0; + dr.g_1 = g_1; + dr.g_2 = g_2; case 3 [g_0, g_1, g_2, g_3] = k_order_perturbation(dr,0,M,options, oo , ['.' ... mexext]); + dr.g_0 = g_0; + dr.g_1 = g_1; + dr.g_2 = g_2; + dr.g_3 = g_3; otherwise error('order > 3 isn''t implemented') end @@ -61,9 +69,3 @@ function [dr,info] = k_order_pert(dr,M,options,oo) dr.ghuu = ghuu; end - if order > 2 - dr.g_0 = g_0; - dr.g_1 = g_1; - dr.g_2 = g_2; - dr.g_3 = g_3; - end \ No newline at end of file diff --git a/matlab/simult_.m b/matlab/simult_.m index 3e4987628..b6a6ce53b 100644 --- a/matlab/simult_.m +++ b/matlab/simult_.m @@ -42,53 +42,69 @@ global M_ options_ it_ y_ = zeros(size(y0,1),iter+M_.maximum_lag); y_(:,1:M_.maximum_lag) = y0; - k1 = [M_.maximum_lag:-1:1]; - k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]); - k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*M_.endo_nbr; - k3 = M_.lead_lag_incidence(1:M_.maximum_lag,:)'; - k3 = find(k3(:)); - k4 = dr.kstate(find(dr.kstate(:,2) < M_.maximum_lag+1),[1 2]); - k4 = k4(:,1)+(M_.maximum_lag+1-k4(:,2))*M_.endo_nbr; - - if iorder == 1 - if ~isempty(dr.ghu) - for i = M_.maximum_lag+1: iter+M_.maximum_lag - tempx1 = y_(dr.order_var,k1); - tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag); - tempx = tempx2(k2); + + if options_.k_order_solver + options_.seed = 77; + ex_ = [zeros(1,M_.exo_nbr); ex_]; + switch options_.order + case 1 + y_ = dynare_simul_(1,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ... + y_(dr.order_var,1),ex_',M_.Sigma_e,options_.seed,dr.ys(dr.order_var),... + zeros(M_.endo_nbr,1),dr.g_1); + case 2 + y_ = dynare_simul_(2,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ... + y_(dr.order_var,1),ex_',M_.Sigma_e,options_.seed,dr.ys(dr.order_var),dr.g_0, ... + dr.g_1,dr.g_2); + case 3 + y_ = dynare_simul_(3,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ... + y_(dr.order_var,1),ex_',M_.Sigma_e,options_.seed,dr.ys(dr.order_var),dr.g_0, ... + dr.g_1,dr.g_2,dr.g_3); + otherwise + error(['order = ' int2str(order) ' isn''t supported']) + end + y_(dr.order_var,:) = y_; + else + k1 = [M_.maximum_lag:-1:1]; + k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]); + k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*M_.endo_nbr; + k3 = M_.lead_lag_incidence(1:M_.maximum_lag,:)'; + k3 = find(k3(:)); + k4 = dr.kstate(find(dr.kstate(:,2) < M_.maximum_lag+1),[1 2]); + k4 = k4(:,1)+(M_.maximum_lag+1-k4(:,2))*M_.endo_nbr; + + if iorder == 1 + if ~isempty(dr.ghu) + for i = M_.maximum_lag+1: iter+M_.maximum_lag + tempx1 = y_(dr.order_var,k1); + tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag); + tempx = tempx2(k2); y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx+dr.ghu* ... ex_(i-M_.maximum_lag,:)'; - k1 = k1+1; + k1 = k1+1; + end + else + for i = M_.maximum_lag+1: iter+M_.maximum_lag + tempx1 = y_(dr.order_var,k1); + tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag); + tempx = tempx2(k2); + y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx; + k1 = k1+1; + end end - else + elseif iorder == 2 for i = M_.maximum_lag+1: iter+M_.maximum_lag tempx1 = y_(dr.order_var,k1); tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag); tempx = tempx2(k2); - y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx; + tempu = ex_(i-M_.maximum_lag,:)'; + tempuu = kron(tempu,tempu); + tempxx = kron(tempx,tempx); + tempxu = kron(tempx,tempu); + y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghs2/2+dr.ghx*tempx+ ... + dr.ghu*tempu+0.5*(dr.ghxx*tempxx+dr.ghuu*tempuu)+dr.ghxu*tempxu; k1 = k1+1; end end - elseif iorder == 2 - for i = M_.maximum_lag+1: iter+M_.maximum_lag - tempx1 = y_(dr.order_var,k1); - tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,M_.maximum_lag); - tempx = tempx2(k2); - tempu = ex_(i-M_.maximum_lag,:)'; - tempuu = kron(tempu,tempu); - tempxx = kron(tempx,tempx); - tempxu = kron(tempx,tempu); - y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghs2/2+dr.ghx*tempx+ ... - dr.ghu*tempu+0.5*(dr.ghxx*tempxx+dr.ghuu*tempuu)+dr.ghxu*tempxu; - k1 = k1+1; - end - elseif iorder == 3 - options_.seed = 77; - ex_ = [zeros(1,M_.exo_nbr); ex_]; - y_ = dynare_simul_(3,dr.nstatic,dr.npred-dr.nboth,dr.nboth,dr.nfwrd,M_.exo_nbr, ... - y_(dr.order_var,1),ex_',M_.Sigma_e,options_.seed,dr.ys(dr.order_var),dr.g_0, ... - dr.g_1,dr.g_2,dr.g_3); - y_(dr.order_var,:) = y_; end % MJ 08/30/02 corrected bug at order 2 \ No newline at end of file diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index b76306f96..6e068a2b1 100644 --- a/matlab/stoch_simul.m +++ b/matlab/stoch_simul.m @@ -27,7 +27,7 @@ function info=stoch_simul(var_list) options_.replic = 1; elseif options_.order == 3 options_.simul = 1; - options_.use_k_order = 1; + options_.k_order_solver = 1; end diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index 678e3ce4f..6d1727f15 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -97,7 +97,7 @@ class ParsingDriver; %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT %token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS %token FLOAT_NUMBER -%token FORECAST +%token FORECAST K_ORDER_SOLVER %token GAMMA_PDF GRAPH CONDITIONAL_VARIANCE_DECOMPOSITION %token HISTVAL HOMOTOPY_SETUP HOMOTOPY_MODE HOMOTOPY_STEPS HP_FILTER HP_NGRID %token IDENTIFICATION INF_CONSTANT INITVAL INITVAL_FILE @@ -804,6 +804,7 @@ stoch_simul_options : o_dr_algo | o_aim_solver | o_partial_information | o_conditional_variance_decomposition + | o_k_order_solver ; symbol_list : symbol_list symbol @@ -1858,6 +1859,7 @@ o_draws_nbr_burn_in_2 : DRAWS_NBR_BURN_IN_2 EQUAL INT_NUMBER {driver.option_num( o_draws_nbr_mean_var_estimate : DRAWS_NBR_MEAN_VAR_ESTIMATE EQUAL INT_NUMBER {driver.option_num("ms.draws_nbr_mean_var_estimate",$3); }; o_draws_nbr_modified_harmonic_mean : DRAWS_NBR_MODIFIED_HARMONIC_MEAN EQUAL INT_NUMBER {driver.option_num("ms.draws_nbr_modified_harmonic_mean",$3); }; o_dirichlet_scale : DIRICHLET_SCALE EQUAL INT_NUMBER {driver.option_num("ms.dirichlet_scale",$3); }; +o_k_order_solver : K_ORDER_SOLVER {driver.option_num("k_order_solver","1"); }; o_chain : CHAIN EQUAL INT_NUMBER { ;}; o_state : STATE EQUAL INT_NUMBER { ;}; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index ba0560ee7..28d2d67cc 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -343,6 +343,7 @@ int sigma_e = 0; posterior_mode {return token::POSTERIOR_MODE; } posterior_mean {return token::POSTERIOR_MEAN; } posterior_median {return token::POSTERIOR_MEDIAN; } +k_order_solver {return token::K_ORDER_SOLVER; } [\$][^$]*[\$] { diff --git a/tests/k_order_perturbation/fs2000k2.mod b/tests/k_order_perturbation/fs2000k2.mod index 07eefba8f..324e87417 100644 --- a/tests/k_order_perturbation/fs2000k2.mod +++ b/tests/k_order_perturbation/fs2000k2.mod @@ -55,7 +55,7 @@ end; steady; -stoch_simul(order=2,use_k_order,periods=1000); +stoch_simul(order=2,k_order_solver,periods=1000); oo1 = load('fs2000k2a_results','oo_'); diff --git a/tests/k_order_perturbation/fs2000k_1.mod b/tests/k_order_perturbation/fs2000k_1.mod index c43f262f2..a611c2eaa 100644 --- a/tests/k_order_perturbation/fs2000k_1.mod +++ b/tests/k_order_perturbation/fs2000k_1.mod @@ -1,4 +1,4 @@ -var m P c e W R k d n l gy_obs gp_obs y dA ; +var m m_1 P P_1 c e W R k d n l gy_obs gp_obs y dA AUXv; varexo e_a e_m; parameters alp bet gam mst rho psi del; @@ -11,10 +11,10 @@ rho = 0.7; psi = 0.787; del = 0.02; -model (use_dll); +model(use_dll); dA = exp(gam+e_a); -log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; --P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +log(m) = (1-rho)*log(mst) + rho*log(m_1(-1))+e_m; +-P/(c(+1)*P(+1)*m)+AUXv(+1)=0; W = l/n; -(psi/(1-psi))*(c*P/(1-n))+l/n = 0; R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; @@ -25,12 +25,17 @@ m-1+d = l; e = exp(e_a); y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); gy_obs = dA*y/y(-1); -gp_obs = (P/P(-1))*m(-1)/dA; +gp_obs = (P/P_1(-1))*m_1(-1)/dA; +m_1 = m; +P_1 = P; +AUXv = bet*P*(alp*exp(-alp*(gam+log(e)))*k(-1)^(alp-1)*n^(1-alp)+(1-del)*exp(-(gam+log(e))))/(c(+1)*P(+1)*m); end; initval; m = mst; +m_1=mst; P = 2.25; +P_1 = 2.25; c = 0.45; e = 1; W = 4; @@ -43,6 +48,7 @@ y = 0.6; gy_obs = exp(gam); gp_obs = exp(-gam); dA = exp(gam); +AUXv = 1; end; shocks; @@ -52,7 +58,7 @@ end; steady; -stoch_simul(order=2,use_k_order,irf=0); +stoch_simul(order=2,k_order_solver,irf=0); if ~exist('fs2000k2_results.mat','file'); error('fs2000k2 must be run first'); @@ -63,27 +69,27 @@ oo1 = load('fs2000k2_results','oo_'); dr0 = oo1.oo_.dr; dr = oo_.dr; -ikr = [15 1:3 16 4:7 12 10 11 13 14]; +ikr = [2:10 1 13:17]; ikc = [1 3 4 2]; ikc2 = [1 3 4 2 9 11 12 10 13 15 16 14 5 7 8 6]; ikc2u = [1 2 5 6 7 8 3 4]; -if max(max(abs(dr0.ghx(ikr,ikc) - dr.ghx(1:end-1,:)))) > 1e-12; +if max(max(abs(dr0.ghx - dr.ghx(ikr,ikc)))) > 1e-12; disp('error in ghx'); end; -if max(max(abs(dr0.ghu(ikr,:) - dr.ghu(1:end-1,:)))) > 1e-12; +if max(max(abs(dr0.ghu - dr.ghu(ikr,:)))) > 1e-12; disp('error in ghu'); end; -if max(max(abs(dr0.ghxx(ikr,ikc2) - dr.ghxx(1:end-1,:)))) > 1e-12; +if max(max(abs(dr0.ghxx - dr.ghxx(ikr,ikc2)))) > 1e-12; disp('error in ghxx'); end; -if max(max(abs(dr0.ghuu(ikr,:) - dr.ghuu(1:end-1,:)))) > 1e-12; +if max(max(abs(dr0.ghuu - dr.ghuu(ikr,:)))) > 1e-12; disp('error in ghuu'); end; -if max(max(abs(dr0.ghxu(ikr,ikc2u) - dr.ghxu(1:end-1,:)))) > 1e-12; +if max(max(abs(dr0.ghxu - dr.ghxu(ikr,ikc2u)))) > 1e-12; disp('error in ghxu'); end; -if max(max(abs(dr0.ghs2(ikr,:) - dr.ghs2(1:end-1,:)))) > 1e-12; +if max(max(abs(dr0.ghs2 - dr.ghs2(ikr,:)))) > 1e-12; disp('error in ghs2'); end;