Commit a615f025 authored by Johannes Pfeifer's avatar Johannes Pfeifer
Browse files
parents 01b51a79 5657dc82
......@@ -145,6 +145,10 @@ public:@;
void walkStochSteady();
TwoDMatrix* calcYCov() const;
const FGSContainer* get_rule_ders() const
{@+ return rule_ders;@+}
const FGSContainer* get_rule_ders_ss() const
{@+ return rule_ders;@+}
protected:@;
void approxAtSteady();
void calcStochShift(Vector& out, double at_sigma) const;
......
......@@ -300,7 +300,7 @@ TwoDMatrix* simulate(emethod em, int np, const Vector& ystart,
@<initialize vectors and subvectors for simulation@>;
@<perform the first step of simulation@>;
@<perform all other steps of simulations@>;
@<add the steady state to all columns of |res|@>;
@<add the steady state to columns of |res|@>;
return res;
}
......@@ -327,10 +327,11 @@ evaluate the polynomial.
eval(em, out, dyu);
@ Also clear. If the result at some period is not finite, we pad the
rest of the matrix with zeros and return immediatelly.
rest of the matrix with zeros.
@<perform all other steps of simulations@>=
for (int i = 1; i < np; i++) {
int i=1;
while (i < np) {
ConstVector ym(*res, i-1);
ConstVector dym(ym, ypart.nstat, ypart.nys());
dy = dym;
......@@ -342,14 +343,16 @@ rest of the matrix with zeros and return immediatelly.
TwoDMatrix rest(*res, i+1, np-i-1);
rest.zeros();
}
return res;
break;
}
i++;
}
@ Even clearer.
@<add the steady state to all columns of |res|@>=
for (int i = 0; i < res->ncols(); i++) {
Vector col(*res, i);
@ Even clearer. We add the steady state to the numbers computed above and leave the padded columns to zero.
@<add the steady state to columns of |res|@>=
for (int j = 0; j < i; j++) {
Vector col(*res, j);
col.add(1.0, ysteady);
}
......
......@@ -63,38 +63,49 @@ function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch)
max_it = 300;
it = 0;
info = 0;
crit = 1+cvg_tol;
A_0 = A1;
X = [];
crit = Inf;
A0_0 = A0;
A1_0 = A1;
A2_0 = A2;
Ahat1 = A1;
if (nargin == 5 && ~isempty(ch) )
A1_0 = A1;
A2_0 = A2;
end
n = length(A0);
id = 1:n;
id0 = 1:n;
id2 = id0+n;
while crit > cvg_tol && it < max_it;
tmp = [A2_0; A0_0]/A1_0;
TMP = tmp*A0_0;
A2_1 = - tmp(id,:)* A2_0;
A_1 = A_0 - TMP(id,:);
A1_1 = A1_0 - tmp(n+id,:) * A2_0 - TMP(id,:);
crit = sum(sum(abs(A0_0)));
A_0 = A_1;
A0_0 = -TMP(n+id,:);
A1_0 = A1_1;
A2_0 = A2_1;
cont = 1;
while cont
tmp = ([A0; A2]/A1)*[A0 A2];
A1 = A1 - tmp(id0,id2) - tmp(id2,id0);
A0 = -tmp(id0,id0);
A2 = -tmp(id2,id2);
Ahat1 = Ahat1 -tmp(id2,id0);
crit = norm(A0,1);
if crit < cvg_tol
% keep iterating until condition on A2 is met
if norm(A2,1) < cvg_tol
cont = 0;
end
elseif isnan(crit) || it == max_it
if crit < cvg_tol
info(1) = 4;
info(2) = log(norm(A2,1));
else
info(1) = 3;
info(2) = log(norm(A1,1));
end
return
end
it = it + 1;
end
if it==max_it
disp(['convergence not achieved after ' int2str(it) ' iterations']);
info = 1;
end
X = -A_0\A0;
X = -Ahat1\A0_0;
if (nargin == 5 && ~isempty(ch) )
%check the solution
res = A0 + A1 * X + A2 * X * X;
res = A0_0 + A1_0 * X + A2_0 * X * X;
if (sum(sum(abs(res))) > cvg_tol)
disp(['the norm residual of the residu ' num2str(res) ' compare to the tolerance criterion ' num2str(cvg_tol)]);
end
......
......@@ -175,31 +175,37 @@ B(:,cols_b) = aa(:,index_c); % Jacobian matrix for contemporaneous endogeneous
C = aa(:,index_p); % Jacobain matrix for led endogeneous variables
info = 0;
info1 = 1;
if task ~= 1 && (DynareOptions.dr_cycle_reduction || DynareOptions.dr_logarithmic_reduction)
if n_current < DynareModel.endo_nbr
if DynareOptions.dr_cycle_reduction
error(['The cycle reduction algorithme can''t be used when the ' ...
'coefficient matrix for current variables is singular'])
'coefficient matrix for current variables isn''t invertible'])
elseif DynareOptions.dr_logarithmic_reduction
error(['The logarithmic reduction algorithme can''t be used when the ' ...
'coefficient matrix for current variables is singular'])
'coefficient matrix for current variables isn''t invertible'])
end
end
end
if DynareOptions.gpu
gpuArray(A1);
gpuArray(B1);
gpuArray(C1);
end
A1 = [aa(row_indx,index_m ) zeros(ndynamic,nfwrd)];
B1 = [aa(row_indx,index_0m) aa(row_indx,index_0p) ];
C1 = [zeros(ndynamic,npred) aa(row_indx,index_p)];
if DynareOptions.dr_cycle_reduction == 1
[ghx, info1] = cycle_reduction(A1, B1, C1, DynareOptions.dr_cycle_reduction_tol);
[ghx, info] = cycle_reduction(A1, B1, C1, DynareOptions.dr_cycle_reduction_tol);
else
[ghx, info1] = logarithmic_reduction(C1, B1, A1, DynareOptions.dr_logarithmic_reduction_tol, DynareOptions.dr_logarithmic_reduction_maxiter);
[ghx, info] = logarithmic_reduction(C1, B1, A1, DynareOptions.dr_logarithmic_reduction_tol, DynareOptions.dr_logarithmic_reduction_maxiter);
end
if info
% cycle_reduction or logarithmic redution failed and set info
return
end
ghx = ghx(:,index_m);
hx = ghx(1:npred+nboth,:);
gx = ghx(1+npred:end,:);
end
if info1 == 1
else
D = zeros(ndynamic+nboth);
E = zeros(ndynamic+nboth);
D(row_indx_de_1,index_d1) = aa(row_indx,index_d);
......@@ -209,11 +215,11 @@ if info1 == 1
E = [-aa(row_indx,[index_m index_0p]) ; [zeros(nboth,nboth+npred) eye(nboth,nboth+nfwrd) ] ];
[err, ss, tt, w, sdim, dr.eigval, info2] = mjdgges(E,D,DynareOptions.qz_criterium);
[err, ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E,D,DynareOptions.qz_criterium);
mexErrCheck('mjdgges', err);
if info2
if info2 == -30
if info1
if info1 == -30
% one eigenvalue is close to 0/0
info(1) = 7;
else
......@@ -255,20 +261,24 @@ if info1 == 1
% derivatives with respect to dynamic state variables
% forward variables
Z = w';
Z11t = Z(indx_stable_root, indx_stable_root)';
Z11 = Z(indx_stable_root, indx_stable_root);
Z21 = Z(indx_explosive_root, indx_stable_root);
Z22 = Z(indx_explosive_root, indx_explosive_root);
if ~isscalar(Z22) && (condest(Z22) > 1e9)
[minus_gx,rc] = linsolve(Z22,Z21);
if rc < 1e-9
% Z22 is near singular
info(1) = 5;
info(2) = condest(Z22);
info(2) = -log(rc);
return;
else
gx = - Z22 \ Z21;
end
gx = -minus_gx;
% predetermined variables
hx = Z11t * inv(tt(indx_stable_root, indx_stable_root)) * ss(indx_stable_root, indx_stable_root) * inv(Z11t);
opts.UT = true;
opts.TRANSA = true;
hx1 = linsolve(tt(indx_stable_root, indx_stable_root),Z11,opts)';
hx2 = linsolve(Z11,ss(indx_stable_root, indx_stable_root)')';
hx = hx1*hx2;
ghx = [hx(k1,:); gx(k2(nboth+1:end),:)];
end
dr.gx = gx;
......
......@@ -93,6 +93,11 @@ if (exist('OCTAVE_VERSION') && ~user_has_octave_forge_package('statistics')) ...
addpath([dynareroot '/missing/nanmean'])
end
% linsolve is missing in Octave
if (exist('OCTAVE_VERSION'))
addpath([dynareroot '/missing/linsolve'])
end
% Add path to MEX files
if exist('OCTAVE_VERSION')
addpath([dynareroot '../mex/octave/']);
......
......@@ -30,7 +30,7 @@ label_width = max(size(deblank(char(headers(1,:),labels)),2))+2;
label_fmt = sprintf('%%-%ds',label_width);
values_length = max(ceil(max(max(log10(abs(values))))),1)+val_precis+1;
values_length = max(ceil(max(max(log10(abs(values(isfinite(values))))))),1)+val_precis+1;
if any(values) < 0
values_length = values_length+1;
end
......
......@@ -521,8 +521,12 @@ options_.graph_save_formats.eps = 1;
options_.graph_save_formats.pdf = 0;
options_.graph_save_formats.fig = 0;
% risky steady state
options_.risky_steadystate = 0;
% use GPU
options_.gpu = 0;
% initialize persistent variables in priordens()
priordens([],[],[],[],[],[],1);
% initialize persistent variables in dyn_first_order_solver()
......
function [x,c] = linsolve(A,B,opts)
% (very imperfect) Clone of Matlab's linsolve.
% Copyright (C) 2010-2011 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
c = [];
x = [];
if nargin == 3
if isfield(opts,'TRANSA')
A = A';
end
end
if nargout == 2
c = rcond(A);
end
x = A\B;
......@@ -78,7 +78,9 @@ chol_S = chol(DynareModel.Sigma_e(i_exo_var,i_exo_var));
for i=1:replic
if ~isempty(DynareModel.Sigma_e)
DynareResults.exo_simul(:,i_exo_var) = randn(DynareOptions.periods,nxs)*chol_S;
% we fill the shocks row wise to have the same values
% independently of the length of the simulation
DynareResults.exo_simul(:,i_exo_var) = randn(nxs,DynareOptions.periods)'*chol_S;
end
y_ = simult_(y0,dr,DynareResults.exo_simul,order);
% elimninating initial value
......
......@@ -95,6 +95,9 @@ void
SimulStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
{
mod_file_struct.simul_present = true;
// The following is necessary to allow shocks+endval+simul in a loop
mod_file_struct.shocks_present_but_simul_not_yet = false;
}
void
......
......@@ -150,9 +150,9 @@ EndValStatement::EndValStatement(const init_values_t &init_values_arg,
void
EndValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
{
if (mod_file_struct.shocks_present)
if (mod_file_struct.shocks_present_but_simul_not_yet)
{
cerr << "ERROR: Putting a \"shocks\" block before an \"endval\" block is not permitted. Please swap the two blocks. This limitation will be removed in the next major release of Dynare." << endl;
cerr << "ERROR: Putting a \"shocks\" block before an \"endval\" block is not permitted. Please swap the two blocks. This limitation will be removed in a future release of Dynare." << endl;
exit(EXIT_FAILURE);
}
}
......
......@@ -196,7 +196,7 @@ void
ShocksStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
{
// Workaround for trac ticket #35
mod_file_struct.shocks_present = true;
mod_file_struct.shocks_present_but_simul_not_yet = true;
// Determine if there is a calibrated measurement error
for (var_and_std_shocks_t::const_iterator it = var_shocks.begin();
......@@ -245,7 +245,7 @@ void
MShocksStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
{
// Workaround for trac ticket #35
mod_file_struct.shocks_present = true;
mod_file_struct.shocks_present_but_simul_not_yet = true;
}
ConditionalForecastPathsStatement::ConditionalForecastPathsStatement(const AbstractShocksStatement::det_shocks_t &paths_arg, const SymbolTable &symbol_table_arg) :
......
......@@ -37,7 +37,7 @@ ModFileStructure::ModFileStructure() :
identification_present(false),
estimation_analytic_derivation(false),
partial_information(false),
shocks_present(false),
shocks_present_but_simul_not_yet(false),
histval_present(false),
k_order_solver(false),
calibrated_measurement_errors(false),
......
......@@ -69,9 +69,11 @@ public:
bool estimation_analytic_derivation;
//! Whether the option partial_information is given to stoch_simul/estimation/osr/ramsey_policy
bool partial_information;
//! Whether a shocks or mshocks block is present
/*! Used for the workaround for trac ticket #35 */
bool shocks_present;
//! Whether a shocks or mshocks block has been parsed and no simul command yet run
/*! Used for the workaround for trac ticket #35. When a simul command is
seen, this flag is cleared in order to allow a sequence
shocks+endval+simul in a loop */
bool shocks_present_but_simul_not_yet;
//! Whether a histval bloc is present
/*! Used for the workaround for trac ticket #157 */
bool histval_present;
......
// See fs2000.mod in the examples/ directory for details on the model
@#define countries = 1:100
var
@#for c in countries
m_@{c} P_@{c} c_@{c} e_@{c} W_@{c} R_@{c} k_@{c} d_@{c} n_@{c} l_@{c} gy_obs_@{c} gp_obs_@{c} y_@{c} dA_@{c}
@#endfor
;
varexo e_a e_m;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
@#for c in countries
dA_@{c} = exp(gam+e_a);
log(m_@{c}) = (1-rho)*log(mst) + rho*log(m_@{c}(-1))+e_m;
-P_@{c}/(c_@{c}(+1)*P_@{c}(+1)*m_@{c})+bet*P_@{c}(+1)*(alp*exp(-alp*(gam+log(e_@{c}(+1))))*k_@{c}^(alp-1)*n_@{c}(+1)^(1-alp)+(1-del)*exp(-(gam+log(e_@{c}(+1)))))/(c_@{c}(+2)*P_@{c}(+2)*m_@{c}(+1))=0;
W_@{c} = l_@{c}/n_@{c};
-(psi/(1-psi))*(c_@{c}*P_@{c}/(1-n_@{c}))+l_@{c}/n_@{c} = 0;
R_@{c} = P_@{c}*(1-alp)*exp(-alp*(gam+e_a))*k_@{c}(-1)^alp*n_@{c}^(-alp)/W_@{c};
1/(c_@{c}*P_@{c})-bet*P_@{c}*(1-alp)*exp(-alp*(gam+e_a))*k_@{c}(-1)^alp*n_@{c}^(1-alp)/(m_@{c}*l_@{c}*c_@{c}(+1)*P_@{c}(+1)) = 0;
c_@{c}+k_@{c} = exp(-alp*(gam+e_a))*k_@{c}(-1)^alp*n_@{c}^(1-alp)+(1-del)*exp(-(gam+e_a))*k_@{c}(-1);
P_@{c}*c_@{c} = m_@{c};
m_@{c}-1+d_@{c} = l_@{c};
e_@{c} = exp(e_a);
y_@{c} = k_@{c}(-1)^alp*n_@{c}^(1-alp)*exp(-alp*(gam+e_a));
gy_obs_@{c} = dA_@{c}*y_@{c}/y_@{c}(-1);
gp_obs_@{c} = (P_@{c}/P_@{c}(-1))*m_@{c}(-1)/dA_@{c};
@#endfor
end;
initval;
@#for c in countries
k_@{c} = 6;
m_@{c} = mst;
P_@{c} = 2.25;
c_@{c} = 0.45;
e_@{c} = 1;
W_@{c} = 4;
R_@{c} = 1.02;
d_@{c} = 0.85;
n_@{c} = 0.19;
l_@{c} = 0.86;
y_@{c} = 0.6;
gy_obs_@{c} = exp(gam);
gp_obs_@{c} = exp(-gam);
dA_@{c} = exp(gam);
@#endfor
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady;
tic;
check;
disp(toc);
tic;
stoch_simul(order=1,dr=cycle_reduction,irf=0,nomoments,noprint);
disp(toc);
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment