v4: global dr_ -> oo_.dr; corrected bug with dr_algo=1

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@608 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
michel 2006-01-18 16:50:33 +00:00
parent 707e1f67bf
commit 6e23ddbf6d
9 changed files with 56 additions and 69 deletions

View File

@ -7,5 +7,5 @@ function ghs2=dr2(ys,dr)
dr.ys = ys;
fh = str2func([M_.fname '_static']);
dr.fbias = 2*feval(fh,dr.ys);
dr=dr1(2,dr,0);
dr=dr1(dr,0);
ghs2 = dr.ghs2;

View File

@ -649,7 +649,7 @@ if (any(bayestopt_.pshape >0 ) & options_.mh_replic) | (any(bayestopt_.pshape >
if options_.mh_replic
metropolis(xparam1,invhess,gend,data,rawdata,bounds);
end
if ~options_.nodiagnostic & options_.mh_replic > 1000 & options_.nblocks > 1
if ~options_.nodiagnostic & options_.mh_replic > 1000 & options_.mh_nblck > 1
McMCDiagnostics;
end
%% Here i discard first half of the draws:

Binary file not shown.

View File

@ -1,7 +1,7 @@
function [A,B,ys,info] = dynare_resolve()
global oo_ dr_
global oo_
[dr_,info] = resol(oo_.steady_state,0);
[oo_.dr,info] = resol(oo_.steady_state,0);
if info(1) > 0
A = [];
@ -10,5 +10,5 @@ function [A,B,ys,info] = dynare_resolve()
return
end
[A,B] = kalman_transition_matrix(dr_);
ys = dr_.ys;
[A,B] = kalman_transition_matrix(oo_.dr);
ys = oo_.dr.ys;

View File

@ -1,6 +1,6 @@
% Copyright (C) 2001 Michel Juillard
%
function [x,cheik] = dynare_solve(func,x,varargin)
function [x,cheik] = dynare_solve(func,x,jacobian_flag,varargin)
global options_
options_ = set_default_option(options_,'solve_algo',2);
@ -11,12 +11,16 @@ function [x,cheik] = dynare_solve(func,x,varargin)
options.MaxFunEvals = 20000;
options.TolFun=1e-8;
options.Display = 'off';
options.Jacobian = 'on';
if jacobian_flag
options.Jacobian = 'on';
else
options.Jacobian = 'on';
end
[x,fval,exitval,output] = fsolve(func,x,options,varargin{:});
if exitval > 0
cheik = 0;
cheik = 0;
else
cheik = 1;
cheik = 1;
end
return
else
@ -26,14 +30,17 @@ function [x,cheik] = dynare_solve(func,x,varargin)
if options_.solve_algo == 1
nn = size(x,1) ;
[x,cheik]=solve1(func,x,1:nn,1:nn,varargin{:});
[x,cheik]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
elseif options_.solve_algo == 2
nn = size(x,1) ;
tolf = eps^(2/3) ;
%fjac = zeros(nn,nn) ;
[fvec,fjac] = feval(func,x,varargin{:});
if jacobian_flag
[fvec,fjac] = feval(func,x,varargin{:});
else
fvec = feval(func,x,varargin{:});
fjac = zeros(nn,nn) ;
end
i = find(~isfinite(fvec));
@ -50,17 +57,20 @@ function [x,cheik] = dynare_solve(func,x,varargin)
return ;
end
%dh = max(abs(x),options_.gstep*ones(nn,1))*eps^(1/3);
%for j = 1:nn
% xdh = x ;
% xdh(j) = xdh(j)+dh(j) ;
% fjac(:,j) = (feval(func,xdh,varargin{:}) - fvec)./dh(j) ;
%end
if ~jacobian_flag
fjac = zeros(nn,nn) ;
dh = max(abs(x),options_.gstep*ones(nn,1))*eps^(1/3);
for j = 1:nn
xdh = x ;
xdh(j) = xdh(j)+dh(j) ;
fjac(:,j) = (feval(func,xdh,varargin{:}) - fvec)./dh(j) ;
end
end
[j1,j2,r,s] = dmperm(fjac);
for i=length(r)-1:-1:1
[x,cheik]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),varargin{:});
[x,cheik]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag,varargin{:});
if cheik
error(sprintf('Solve block = %d check = %d\n',i,cheik));
end

View File

@ -57,7 +57,7 @@ if info(1)
end
if options_.dr_algo == 1 & options_.order > 1
dr.ys = dynare_solve('dr2',ys,dr);
dr.ys = dynare_solve('dr2',ys,0,dr);
dr.fbias = 2*feval([M_.fname '_static'],dr.ys,oo_.exo_steady_state);
[dr, info1] = dr1(dr,check_flag);
if info1(1)

View File

@ -1,6 +1,6 @@
% Copyright (C) 2001 Michel Juillard
%
function [x,check] = solve1(func,x,j1,j2,varargin)
function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin)
global M_ options_ fjac
@ -37,17 +37,19 @@ function [x,check] = solve1(func,x,j1,j2,varargin)
stpmax = stpmx*max([sqrt(x'*x);nn]) ;
first_time = 1;
for its = 1:maxit
%dh = max(abs(x(j2)),options_.gstep*ones(nn,1))*eps^(1/3);
%for j = 1:nn
% xdh = x ;
% xdh(j2(j)) = xdh(j2(j))+dh(j) ;
% t = feval(func,xdh,varargin{:});
% fjac(:,j) = (t(j1) - fvec)./dh(j) ;
% g(j) = fvec'*fjac(:,j) ;
%end
[fvec,fjac] = feval(func,x,varargin{:});
if jacobian_flag
[fvec,fjac] = feval(func,x,varargin{:});
else
dh = max(abs(x(j2)),options_.gstep*ones(nn,1))*eps^(1/3);
for j = 1:nn
xdh = x ;
xdh(j2(j)) = xdh(j2(j))+dh(j) ;
t = feval(func,xdh,varargin{:});
fjac(:,j) = (t(j1) - fvec)./dh(j) ;
g(j) = fvec'*fjac(:,j) ;
end
end
fvec = fvec(j1);
fjac = fjac(j1,j2);
g = (fvec'*fjac)';

View File

@ -18,7 +18,7 @@ function steady_()
if exist([M_.fname '_steadystate'])
[oo_.steady_state,check] = feval([M_.fname '_steadystate'],oo_.steady_state,x);
else
[oo_.staedy_state,check] = dynare_solve([M_.fname '_static'],oo_.steady_state,x);
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],oo_.steady_state,1,x);
end
if check ~= 0

View File

@ -1,7 +1,7 @@
% Copyright (C) 2001 Michel Juillard
%
function info=stoch_simul(var_list)
global M_ options_ oo_ dr_
global M_ options_ oo_
global it_
options_ = set_default_option(options_,'TeX',0);
@ -42,7 +42,7 @@ global it_
check_model;
[dr_, info] = resol(oo_.steady_state,0);
[oo_.dr, info] = resol(oo_.steady_state,0);
if info(1)
print_info(info);
@ -56,21 +56,21 @@ global it_
disp([' Number of variables: ' int2str(M_.endo_nbr)])
disp([' Number of stochastic shocks: ' int2str(M_.exo_nbr)])
disp([' Number of state variables: ' ...
int2str(length(find(dr_.kstate(:,2) <= M_.maximum_lag+1)))])
int2str(length(find(oo_.dr.kstate(:,2) <= M_.maximum_lag+1)))])
disp([' Number of jumpers: ' ...
int2str(length(find(dr_.kstate(:,2) == M_.maximum_lag+2)))])
disp([' Number of static variables: ' int2str(dr_.nstatic)])
int2str(length(find(oo_.dr.kstate(:,2) == M_.maximum_lag+2)))])
disp([' Number of static variables: ' int2str(oo_.dr.nstatic)])
my_title='MATRIX OF COVARIANCE OF EXOGENOUS SHOCKS';
labels = deblank(M_.exo_names);
headers = strvcat('Variables',labels);
lh = size(labels,2)+2;
table(my_title,headers,labels,M_.Sigma_e,lh,10,6);
disp(' ')
disp_dr(dr_,options_.order,var_list);
disp_dr(oo_.dr,options_.order,var_list);
end
if options_.simul == 0 & options_.nomoments == 0
disp_th_moments(dr_,var_list);
disp_th_moments(oo_.dr,var_list);
elseif options_.simul == 1
if options_.periods == 0
error('STOCH_SIMUL error: number of periods for the simulation isn''t specified')
@ -80,7 +80,7 @@ global it_
' than the number of observations to be DROPed'])
return
end
oo_.endo_simul = simult(repmat(dr_.ys,1,M_.maximum_lag),dr_);
oo_.endo_simul = simult(repmat(oo_.dr.ys,1,M_.maximum_lag),oo_.dr);
dyn2vec;
if options_.nomoments == 0
disp_moments(oo_.endo_simul,var_list);
@ -130,7 +130,7 @@ global it_
end
for i=1:M_.exo_nbr
if SS(i,i) > 1e-13
y=irf(dr_,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ...
y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ...
options_.replic, options_.order);
if options_.relative_irf
y = 100*y/cs(i,i);
@ -273,28 +273,3 @@ global it_
fclose(fidTeX);
end
end
% 02/20/01 MJ oo_.steady_state removed from calling sequence for simult (all in dr_)
% 02/23/01 MJ added dyn2vec()
% 06/24/01 MJ steady -> steadoo_.endo_simul
% 08/28/02 MJ added var_list
% 10/09/02 MJ no simulation and theoretical moments for order 1
% 10/14/02 MJ added plot of IRFs
% 10/30/02 MJ options_ are now a structure
% 01/01/03 MJ added dr_algo
% 01/09/03 MJ set default values for options_ (correct absence of autocorr
% when order == 1)
% 01/12/03 MJ removed call to steadoo_.endo_simul as already checked in resol()
% 02/09/03 MJ oo_.steady_state reset with value declared in initval after computations
% 02/18/03 MJ removed above change. oo_.steady_state shouldn't be affected by
% computations in this function
% new option SIMUL computes a stochastic simulation and save
% results in oo_.endo_simul and via dyn2vec
% 04/03/03 MJ corrected bug for simulation with M_.maximum_lag > 1
% 05/20/03 MJ eliminates exogenous shocks with 0 variance
% 05/20/03 MJ don't plot IRF if variation < 1e-10
% 11/14/03 MJ corrected bug on number of replications for IRF when
% order=2
% 11/22/03 MJ replaced IRFs by orthogonalized IRFs
% 08/30/04 SA The maximum number of plots is not constrained for the IRFs and
% all the plots are saved in *.eps, *.pdf and *.fig files (added
% 09/03/04 SA Tex output for IRFs added