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

Merge pull request #1200 from rattoma/gsa 

identification and slice
time-shift
Houtan Bastani 2016-05-23 11:59:26 +02:00
parent 389f43019d 3e8dbcfc7e
commit 25f6d81cb7
7 changed files with 63 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
matlab/check_posterior_sampler_options.m
View File

@ -230,10 +230,10 @@ if init,
end
case 'mode_files'
% for multimodal posteriors provide a list of mode files,
% one per mode. With this info, the code will automatically
% set the <mode> option. The mode files need only to
% contain the xparam1 variable.
% for multimodal posteriors provide the name of
% a file containing a variable array xparams = [nparam * nmodes]
% one column per mode. With this info, the code will automatically
% set the <mode> option.
% This will automatically trigger <rotated>
% default = []
posterior_sampler_options.mode_files = options_list{i,2};
@ -326,9 +326,12 @@ if init,
if ~isempty(posterior_sampler_options.mode_files), % multimodal case
modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains
for j=1:length( modes ),
load(modes{j}, 'xparam1')
mode(j).m=xparam1;
load(modes, 'xparams')
if size(xparams,2)<2,
error(['check_posterior_sampler_options:: Variable xparams loaded in file <' modes '> has size [' int2str(size(xparams,1)) 'x' int2str(size(xparams,2)) ']: it must contain at least two columns, to allow multi-modal sampling.'])
end
for j=1:size(xparams,2),
mode(j).m=xparams(:,j);
end
posterior_sampler_options.mode = mode;
posterior_sampler_options.rotated = 1;

4
matlab/dsge_likelihood.m
View File

@ -546,10 +546,10 @@ if analytic_derivation,
end
if full_Hess,
[dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv);
[dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv);
clear dum dum2;
else
[dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv);
[dum, DT, DOm, DYss] = getH(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,iv);
end
else
DT = derivatives_info.DT(iv,iv,:);

16
matlab/dynare_identification.m
View File

@ -187,10 +187,18 @@ if prior_exist,
name = bayestopt_.name;
name_tex = char(M_.exo_names_tex(indexo,:),M_.param_names_tex(indx,:));
offset = estim_params_.nvx;
offset = offset + estim_params_.nvn;
offset = offset + estim_params_.ncx;
offset = offset + estim_params_.ncn;
if estim_params_.nvn || estim_params_.ncn,
error('Identification does not support measurement errors. Instead, define them explicitly in measurement equations in model definition.')
else
offset = estim_params_.nvx;
%offset = offset + estim_params_.nvn;
offset = offset + estim_params_.ncx;
if estim_params_.ncx
options_ident.analytic_derivation=0;
options_ident.analytic_derivation_mode=-1;
end
%offset = offset + estim_params_.ncn;
end
else
indx = [1:M_.param_nbr];
indexo = [1:M_.exo_nbr];

18
matlab/getH.m
View File

@ -1,5 +1,5 @@
function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,options_,kronflag,indx,indexo,iv)
% function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,options_,kronflag,indx,indexo,iv)
function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, estim_params_,M_,oo_,options_,kronflag,indx,indexo,iv)
% function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, estim_params_,M_,oo_,options_,kronflag,indx,indexo,iv)
% computes derivative of reduced form linear model w.r.t. deep params
%
% Inputs:
@ -49,16 +49,16 @@ function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,options_,kro
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<6 || isempty(kronflag)
if nargin<7 || isempty(kronflag)
kronflag = 0;
end
if nargin<7 || isempty(indx)
if nargin<8 || isempty(indx)
indx = [];
end
if nargin<8 || isempty(indexo)
if nargin<9 || isempty(indexo)
indexo = [];
end
if nargin<9 || isempty(iv)
if nargin<10 || isempty(iv)
iv = (1:length(A))';
end
@ -96,7 +96,7 @@ if kronflag==-1, % perturbation
end
if nargout>5,
H2 = hessian_sparse('thet2tau',[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)], ...
options_.gstep,M_, oo_, indx,indexo,0,[],[],[],iv);
options_.gstep,estim_params_,M_, oo_, indx,indexo,0,[],[],[],iv);
H2ss = zeros(m1,tot_param_nbr,tot_param_nbr);
iax=find(triu(rand(tot_param_nbr,tot_param_nbr)));
H2 = H2(:,iax);
@ -126,7 +126,7 @@ if kronflag==-2,
if nargout>5,
[residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
M_.params, oo_.dr.ys, 1);
g22 = hessian_sparse('thet2tau',[M_.params(indx)],options_.gstep,M_, oo_, indx,[],-1);
g22 = hessian_sparse('thet2tau',[M_.params(indx)],options_.gstep,estim_params_,M_, oo_, indx,[],-1);
H2ss=full(g22(1:M_.endo_nbr,:));
H2ss = reshape(H2ss,[M_.endo_nbr param_nbr param_nbr]);
for j=1:M_.endo_nbr,
@ -147,7 +147,7 @@ if kronflag==-2,
[residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
M_.params, oo_.dr.ys, 1);
end
gp = fjaco('thet2tau',[M_.params(indx)],M_, oo_, indx,[],-1);
gp = fjaco('thet2tau',[M_.params(indx)],estim_params_,M_, oo_, indx,[],-1);
Hss=gp(1:M_.endo_nbr,:);
gp=gp(M_.endo_nbr+1:end,:);
gp = reshape(gp,[size(g1) param_nbr]);

25
matlab/getJJ.m
View File

@ -1,5 +1,5 @@
function [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr)
% function [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr)
function [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, estim_params_, M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr)
% function [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, estim_params_, M_,oo_,options_,kronflag,indx,indexo,mf,nlags,useautocorr)
% computes derivatives of 1st and 2nd order moments of observables with
% respect to estimated parameters
%
@ -54,16 +54,16 @@ function [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo_,options_,kronflag,
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<7 || isempty(indx)
if nargin<8 || isempty(indx)
% indx = [1:M_.param_nbr];
end,
if nargin<8 || isempty(indexo)
if nargin<9 || isempty(indexo)
indexo = [];
end,
if nargin<10 || isempty(nlags)
if nargin<11 || isempty(nlags)
nlags=3;
end
if nargin<11 || isempty(useautocorr)
if nargin<12 || isempty(useautocorr)
useautocorr=0;
end
@ -73,15 +73,16 @@ warning('off','MATLAB:divideByZero')
if kronflag == -1,
fun = 'thet2tau';
params0 = M_.params;
JJ = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],M_, oo_, indx,indexo,1,mf,nlags,useautocorr);
para0 = get_all_parameters(estim_params_, M_);
JJ = fjaco(fun,para0,estim_params_,M_, oo_, indx,indexo,1,mf,nlags,useautocorr);
M_.params = params0;
params0 = M_.params;
H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],M_, oo_, indx,indexo,0,mf,nlags,useautocorr);
H = fjaco(fun,para0,estim_params_,M_, oo_, indx,indexo,0,mf,nlags,useautocorr);
M_.params = params0;
params0 = M_.params;
gp = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],M_, oo_, indx,indexo,-1);
gp = fjaco(fun,para0,estim_params_,M_, oo_, indx,indexo,-1);
M_.params = params0;
offset = length(indexo);
offset = length(para0)-length(indx);
gp = gp(:,offset+1:end);
dYss = H(1:M_.endo_nbr,offset+1:end);
dA = reshape(H(M_.orig_endo_nbr+[1:numel(A)],:),[size(A),size(H,2)]);
@ -92,7 +93,7 @@ if kronflag == -1,
assignin('base','M_', M_);
assignin('base','oo_', oo_);
else
[H, dA, dOm, dYss, gp] = getH(A, B, M_,oo_,options_,kronflag,indx,indexo);
[H, dA, dOm, dYss, gp] = getH(A, B, estim_params_,M_,oo_,options_,kronflag,indx,indexo);
gp = reshape(gp,size(gp,1)*size(gp,2),size(gp,3));
gp = [dYss; gp];
% if isempty(H),
@ -194,4 +195,4 @@ gam = [oo_.dr.ys(oo_.dr.order_var(mf)); gam];
% if useautocorr,
warning('on','MATLAB:divideByZero')
% end
% end

7
matlab/identification_analysis.m
View File

@ -80,7 +80,7 @@ if info(1)==0,
oo_.dr.ys, 1);
vg1 = [oo_.dr.ys(oo_.dr.order_var); vec(g1)];
[JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,kron_flag,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
[JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, estim_params_, M_,oo0,options_,kron_flag,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
derivatives_info.DT=dA;
derivatives_info.DOm=dOm;
derivatives_info.DYss=dYss;
@ -95,7 +95,7 @@ if info(1)==0,
disp('The number of moments with non-zero derivative is smaller than the number of parameters')
disp(['Try increasing ar = ', int2str(nlags+1)])
nlags=nlags+1;
[JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,kron_flag,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
[JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, estim_params_, M_,oo0,options_,kron_flag,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
derivatives_info.DT=dA;
derivatives_info.DOm=dOm;
derivatives_info.DYss=dYss;
@ -127,6 +127,9 @@ if info(1)==0,
else
normaliz1=[];
end
if ~isempty(estim_params_.corrx),
normaliz1 = [normaliz1 estim_params_.corrx(:,8)']; % normalize with prior standard deviation
end
if ~isempty(estim_params_.param_vals),
normaliz1 = [normaliz1 estim_params_.param_vals(:,7)']; % normalize with prior standard deviation
end

20
matlab/thet2tau.m
View File

@ -1,4 +1,4 @@
function tau = thet2tau(params, M_, oo_, indx, indexo, flagmoments,mf,nlags,useautocorr,iv)
function tau = thet2tau(params, estim_params_, M_, oo_, indx, indexo, flagmoments,mf,nlags,useautocorr,iv)
%
% Copyright (C) 2011-2012 Dynare Team
@ -25,20 +25,24 @@ if nargin==1,
indexo = [];
end
if nargin<6,
if nargin<7,
flagmoments=0;
end
if nargin<9 || isempty(useautocorr),
if nargin<10 || isempty(useautocorr),
useautocorr=0;
end
if nargin<10 || isempty(iv),
if nargin<11 || isempty(iv),
iv=[1:M_.endo_nbr];
end
M_.params(indx) = params(length(indexo)+1:end);
if ~isempty(indexo)
M_.Sigma_e(indexo,indexo) = diag(params(1:length(indexo)).^2);
if length(params)>length(indx),
M_ = set_all_parameters(params,estim_params_,M_);
else
M_.params(indx) = params;
end
% if ~isempty(indexo)
% M_.Sigma_e(indexo,indexo) = diag(params(1:length(indexo)).^2);
% end
[A,B,tele,tubbies,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
if flagmoments==0,
ys=oo_.dr.ys(oo_.dr.order_var);
@ -71,4 +75,4 @@ else
tau = [tau;vec(dum(mf,mf))];
end
tau = [ oo_.dr.ys(oo_.dr.order_var(mf)); tau];
end
end