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

v4: bugs in trend_coeff, steady state; changed handling of penalty by csminwel 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@590 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
michel 2006-01-08 08:39:00 +00:00
parent 70dffd1518
commit 71a4fd49c8
7 changed files with 86 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
matlab/DsgeLikelihood.m
View File

@ -1,13 +1,13 @@
function [fval,cost_flag,ys,trend_coeff] = DsgeLikelihood(xparam1,gend,data)
function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data)
% stephane.adjemian@cepremap.cnrs.fr [09-07-2004]
%
% Adapted from mj_optmumlik.m
global bayestopt_ estim_params_ options_ trend_coeff_ M_ oo_
global dr1_test_
global bayestopt_ estim_params_ options_ trend_coeff_ M_ oo_ xparam1_test
fval = [];
ys = [];
trend_coeff = [];
xparam1_test = xparam1;
cost_flag = 1;
nobs = size(options_.varobs,1);
%------------------------------------------------------------------------------
@ -15,13 +15,13 @@ function [fval,cost_flag,ys,trend_coeff] = DsgeLikelihood(xparam1,gend,data)
%------------------------------------------------------------------------------
if options_.mode_compute ~= 1 & any(xparam1 < bayestopt_.lb)
k = find(xparam1 < bayestopt_.lb);
fval = bayestopt_.penalty*min(1e3,exp(sum(bayestopt_.lb(k)-xparam1(k))));
fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
cost_flag = 0;
return;
end
if options_.mode_compute ~= 1 & any(xparam1 > bayestopt_.ub)
k = find(xparam1 > bayestopt_.ub);
fval = bayestopt_.penalty*min(1e3,exp(sum(xparam1(k)-bayestopt_.ub(k))));
fval = bayestopt_.penalty+sum((xparam1(k)-bayestopt_.ub(k)).^2);
cost_flag = 0;
return;
end
@ -50,9 +50,12 @@ function [fval,cost_flag,ys,trend_coeff] = DsgeLikelihood(xparam1,gend,data)
if testQ %% The variance-covariance matrix of the structural innovations is not definite positive.
%% We have to compute the eigenvalues of this matrix in order to build the penalty.
a = diag(eig(Q));
fval = bayestopt_.penalty*min(1e3,exp(sum(-a(a<=0))));
cost_flag = 0;
return
k = find(a < 0);
if k > 0
fval = bayestopt_.penalty+sum(-a(k));
cost_flag = 0;
return
end
end
offset = offset+estim_params_.ncx;
end
@ -66,24 +69,11 @@ function [fval,cost_flag,ys,trend_coeff] = DsgeLikelihood(xparam1,gend,data)
[CholH,testH] = chol(H);
if testH
a = diag(eig(H));
if nobs == estim_params_.nvn
fval = bayestopt_.penalty*min(1e3,exp(sum(-a(a<=0))));
k = find(a < 0);
if k > 0
fval = bayestopt_.penalty+sum(-a(k));
cost_flag = 0;
return
else
if sum(abs(a)<crit) == nobs-estim_params_.nvn
if any(a<0)
fval = bayestopt_.penalty*min(1e3,exp(sum(-a(a<0))));
cost_flag = 0;
return
else
% All is fine, there's nothing to do here...
end
else
fval = bayestopt_.penalty*min(1e3,exp(sum(-a(a<=0))));
cost_flag = 0;
return
end
end
end
offset = offset+estim_params_.ncn;
@ -96,17 +86,13 @@ function [fval,cost_flag,ys,trend_coeff] = DsgeLikelihood(xparam1,gend,data)
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
[T,R,SteadyState] = dynare_resolve;
if dr1_test_(1) == 1
fval = bayestopt_.penalty*min(1e3,exp(dr1_test_(2)));
[T,R,SteadyState,info] = dynare_resolve;
if info(1) == 1 | info(1) == 2 | info(1) == 5
fval = bayestopt_.penalty+1;
cost_flag = 0;
return
elseif dr1_test_(1) == 2
fval = bayestopt_.penalty*min(1e3,exp(dr1_test_(2)));
cost_flag = 0;
return
elseif dr1_test_(1) == 3
fval = bayestopt_.penalty*min(1e3,exp(dr1_test_(2)));
elseif info(1) == 3 | info(1) == 4 | info(1) == 20
fval = bayestopt_.penalty+info(2)^2;
cost_flag = 0;
return
end
@ -120,9 +106,9 @@ function [fval,cost_flag,ys,trend_coeff] = DsgeLikelihood(xparam1,gend,data)
for i=1:nobs
trend_coeff(i) = evalin('base',bayestopt_.trend_coeff{i});
end
trend = constant*ones(1,gend)+trend_coeff*(1:gend);
trend = repmat(constant,1,gend)+trend_coeff*[1:gend];
else
trend = constant*ones(1,gend);
trend = repmat(constant,1,gend);
end
start = options_.presample+1;
np = size(T,1);

6
matlab/csminwel.m
View File

@ -16,7 +16,8 @@ function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,va
% hessian update, respectively. (When the routine hits certain kinds of difficulty, it
% write g2.mat and g3.mat as well. If all were written at about the same time, any of them
% may be a decent starting point. One can also start from the one with best function value.)
[nx,no]=size(x0);
global bayestopt_
[nx,no]=size(x0);
nx=max(nx,no);
Verbose=1;
NumGrad= isempty(grad);
@ -60,7 +61,8 @@ f=f0;
H=H0;
cliff=0;
while ~done
g1=[]; g2=[]; g3=[];
bayestopt_.penalty = f;
g1=[]; g2=[]; g3=[];
%addition fj. 7/6/94 for control
disp('-----------------')
disp('-----------------')

49
matlab/dr1.m
View File

@ -139,6 +139,32 @@ sdyn = M_.endo_nbr - nstatic;
k0 = M_.lead_lag_incidence(M_.maximum_lag+1,order_var);
k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_lag+1),:);
b = jacobia_(:,k0);
if M_.maximum_lead == 0; % backward models
a = jacobia_(:,nonzeros(k1'));
dr.ghx = zeros(size(a));
m = 0;
for i=M_.maximum_lag:-1:1
k = nonzeros(M_.lead_lag_incidence(i,order_var));
dr.ghx(:,m+[1:length(k)]) = -b\a(:,k);
m = m+length(k);
end
if M_.exo_nbr
dr.ghu = -b\fu;
end
dr.eigval = eig(transition_matrix(dr));
dr.rank = 0;
if any(abs(dr.eigval) > options_.qz_criterium)
temp = sort(abs(dr.eigval));
nba = nnz(abs(dr.eigval) > options_.qz_criterium);
temp = temp(nd-nba+1:nd)-1-options_.qz_criterium;
info(1) = 3;
info(2) = temp'*temp;
end
return;
end
%forward--looking models
if nstatic > 0
[Q,R] = qr(b(:,1:nstatic));
aa = Q'*jacobia_;
@ -160,29 +186,6 @@ if M_.exo_nbr
end
clear aa;
if M_.maximum_lead == 0; % backward model
dr.ghx = zeros(size(a));
m = 0;
for i=M_.maximum_lag:-1:1
k = nonzeros(M_.lead_lag_incidence(i,order_var));
dr.ghx(:,m+[1:length(k)]) = -b\a(:,k);
m = m+length(k);
end
if M_.exo_nbr
dr.ghu = -fu;
end
dr.eigval = eig(transition_matrix(dr));
dr.rank = 0;
if any(abs(dr.eigval) > options_.qz_criterium)
temp = sort(abs(dr.eigval));
nba = nnz(abs(dr.eigval) > options_.qz_criterium);
temp = temp(nd-nba+1:nd)-1-options_.qz_criterium;
info(1) = 3;
info(2) = temp'*temp;
end
return;
end
% buildind D and E
d = zeros(nd,nd) ;
e = d ;

12
matlab/dynare_estimation.m
View File

@ -1,7 +1,7 @@
function dynare_estimation(var_list_)
global M_ options_ oo_ estim_params_
global bayestopt_ trend_coeff_
global bayestopt_
options_.varlist = var_list_;
options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1);
@ -43,7 +43,7 @@ options_ = set_default_option(options_,'relative_irf',0);
options_ = set_default_option(options_,'order',1);
options_ = set_default_option(options_,'ar',5);
options_ = set_default_option(options_,'dr_algo',0);
options_ = set_default_option(options_,'linear',1);
options_ = set_default_option(options_,'linear',0);
options_ = set_default_option(options_,'drop',0);
options_ = set_default_option(options_,'replic',1);
options_ = set_default_option(options_,'hp_filter',0);
@ -86,16 +86,18 @@ ub = bounds(:,2);
bayestopt_.lb = lb;
bayestopt_.ub = ub;
if isempty(trend_coeff_)
if isempty(options_.trend_coeffs)
bayestopt_.with_trend = 0;
else
bayestopt_.with_trend = 1;
bayestopt_.trend_coeff = {};
trend_coeffs = options_.trend_coeffs;
nt = length(trend_coeffs);
for i=1:n_varobs
if i > length(trend_coeff_) | isempty(trend_coeff_{i})
if i > length(trend_coeffs) | isempty(trend_coeffs{i})
bayestopt_.trend_coeff{i} = '0';
else
bayestopt_.trend_coeff{i} = trend_coeff_{i};
bayestopt_.trend_coeff{i} = trend_coeffs{i};
end
end
end

BIN
matlab/dynare_m.exe
View File

Binary file not shown.

31
matlab/dynare_resolve.m
View File

@ -1,17 +1,14 @@
function [A,B,ys] = dynare_resolve()
global dr1_test_ oo_
dr = resol1(oo_.steady_state,0,1,1);
oo_.dr = dr;
if dr1_test_(1) > 0
A = [];
B = [];
ys = [];
return
end
[A,B] = kalman_transition_matrix(dr);
ys = dr.ys;
% SA 01-18-2005 Added line five (oo_.dr = dr;)
function [A,B,ys,info] = dynare_resolve()
global oo_ dr_
[dr_,info] = resol(oo_.steady_state,0);
if info(1) > 0
A = [];
B = [];
ys = [];
return
end
[A,B] = kalman_transition_matrix(dr_);
ys = dr_.ys;

20
matlab/steady_.m
View File

@ -12,14 +12,22 @@ function steady_()
temp = oo_.exo_simul ;
oo_.exo_simul = ones(xlen,1)*transpose(oo_.exo_steady_state);
[oo_.steady_state,cheik] = dynare_solve([M_.fname '_static'],x,oo_.exo_simul);
if M_.exo_det_nbr > 0
tempdet = oo_.exo_det_simul ;
oo_.exo_det_simul = repmat(oo_.exo_det_steady_state',M_.maximum_lag+1,1) ;
end
if cheik ~= 0
if exist([M_.fname '_steadystate'])
[oo_.steady_state,check] = feval([M_.fname '_steadystate'],x);
else
[oo_.staedy_state,check] = dynare_solve([M_.fname '_static'],x);
end
if check ~= 0
error('STEADY: convergence problems')
end
if M_.exo_det_nbr > 0
oo_.exo_det_simul = tempdet;
end
oo_.exo_simul = temp ;
% 06/24/01 MJ: steady_ no results printer; steady with printed results
% 07/31/03 MJ: in case of convergence problem steady stops with an error
% 01/22/05 SA: check --> cheik (to avoid confusion with function check.m)