Commit 41155369 authored by Houtan Bastani's avatar Houtan Bastani
Browse files

Merge pull request #955 from JohannesPfeifer/identification_fixes

Various fixes to GSA and Identification
parents aaee660f 249c2ca7
......@@ -5358,6 +5358,7 @@ decomposition of the above k-step ahead filtered values. Stores results in @code
@item diffuse_filter
@anchor{diffuse_filter}
Uses the diffuse Kalman filter (as described in
@cite{Durbin and Koopman (2012)} and @cite{Koopman and Durbin
(2003)}) to estimate models with non-stationary observed variables.
......@@ -7238,8 +7239,8 @@ for identification analysis. Default: @code{0}
@item ar = @var{INTEGER}
Maximum number of lags for moments in identification analysis. Default: @code{1}
@item lik_init = @var{INTEGER}
@xref{lik_init}.
@item diffuse_filter = @var{INTEGER}
@xref{diffuse_filter}.
@end table
......
......@@ -59,8 +59,24 @@ options_ident = set_default_option(options_ident,'periods',300);
options_ident = set_default_option(options_ident,'replic',100);
options_ident = set_default_option(options_ident,'advanced',0);
options_ident = set_default_option(options_ident,'normalize_jacobians',1);
%Deal with non-stationary cases
if isfield(options_ident,'diffuse_filter') %set lik_init and options_
options_ident.lik_init=3; %overwrites any other lik_init set
options_.diffuse_filter=options_ident.diffuse_filter; %make options_ inherit diffuse filter; will overwrite any conflicting lik_init in dynare_estimation_init
else
if options_.diffuse_filter==1 %warning if estimation with diffuse filter was done, but not passed
warning('IDENTIFICATION:: Previously the diffuse_filter option was used, but it was not passed to the identification command. This may result in problems if your model contains unit roots.')
end
if isfield(options_ident,'lik_init')
if options_ident.lik_init==3 %user specified diffuse filter using the lik_init option
options_ident.analytic_derivation=0; %diffuse filter not compatible with analytic derivation
end
end
end
options_ident = set_default_option(options_ident,'lik_init',1);
options_ident = set_default_option(options_ident,'analytic_derivation',1);
if isfield(options_ident,'nograph'),
options_.nograph=options_ident.nograph;
end
......@@ -70,6 +86,9 @@ end
if isfield(options_ident,'graph_format'),
options_.graph_format=options_ident.graph_format;
end
if isfield(options_ident,'prior_trunc'),
options_.prior_trunc=options_ident.prior_trunc;
end
if options_ident.gsa_sample_file,
GSAFolder = checkpath('gsa',M_.dname);
......@@ -122,8 +141,7 @@ options_.prior_mc = options_ident.prior_mc;
options_.options_ident = options_ident;
options_.Schur_vec_tol = 1.e-8;
options_.nomoments=0;
options_.analytic_derivation=1;
options_ = set_default_option(options_,'analytic_derivation',1);
options_ = set_default_option(options_,'datafile','');
options_.mode_compute = 0;
options_.plot_priors = 0;
......
......@@ -52,6 +52,7 @@ else
end
i = 1;
c = zeros(n,1,p);
c1 = zeros(n,1,p);
while i < m
if t(i+1,i) == 0
if i == 1
......
......@@ -19,10 +19,18 @@ 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), kronflag = 0; end
if nargin<7 || isempty(indx), indx = [1:M_.param_nbr]; end,
if nargin<8 || isempty(indexo), indexo = []; end,
if nargin<9 || isempty(iv), iv = (1:length(A))'; end,
if nargin<6 || isempty(kronflag)
kronflag = 0;
end
if nargin<7 || isempty(indx)
indx = [];
end
if nargin<8 || isempty(indexo)
indexo = [];
end
if nargin<9 || isempty(iv)
iv = (1:length(A))';
end
[I,J]=find(M_.lead_lag_incidence');
yy0=oo_.dr.ys(I);
......
......@@ -17,10 +17,18 @@ 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), indx = [1:M_.param_nbr];, end,
if nargin<8 || isempty(indexo), indexo = [];, end,
if nargin<10 || isempty(nlags), nlags=3; end,
if nargin<11 || isempty(useautocorr), useautocorr=0; end,
if nargin<7 || isempty(indx)
% indx = [1:M_.param_nbr];
end,
if nargin<8 || isempty(indexo)
indexo = [];
end,
if nargin<10 || isempty(nlags)
nlags=3;
end
if nargin<11 || isempty(useautocorr)
useautocorr=0;
end
% if useautocorr,
warning('off','MATLAB:divideByZero')
......
......@@ -64,7 +64,7 @@ if opt_gsa.load_ident_files==0,
if opt_gsa.morris==2,
pdraws = dynare_identification(options_.options_ident,[lpmatx lpmat(istable,:)]);
% [pdraws, TAU, GAM] = dynare_identification(options_.options_ident,[lpmatx lpmat(istable,:)]);
if max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0,
if ~isempty(pdraws) && max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0,
disp(['Sample check OK ', num2str(max(max(abs(pdraws-[lpmatx lpmat(istable,:)]))))]),
clear pdraws;
end
......
......@@ -20,7 +20,7 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr)
global options_ M_ estim_params_ oo_
[nr1, nc1, nsam] = size(mm);
nobs=size(options_.varobs,1);
nobs=size(options_.varobs,2);
disp('Computing theoretical moments ...')
h = dyn_waitbar(0,'Theoretical moments ...');
vdec = zeros(nobs,M_.exo_nbr,nsam);
......@@ -37,7 +37,7 @@ global options_ M_ estim_params_ oo_
cc(:,:,j)=triu(corr);
dum=[];
for i=1:options_.ar
dum=[dum, autocorr{i}];
dum=[dum, autocorr{i}];
end
ac(:,:,j)=dum;
dyn_waitbar(j/nsam,h)
......
......@@ -42,7 +42,7 @@ function pdraw = prior_draw_gsa(init,rdraw)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global bayestopt_ options_
global bayestopt_ options_ estim_params_ M_
persistent npar pshape p6 p7 p3 p4 lbcum ubcum
if init
......@@ -55,7 +55,18 @@ if init
pdraw = zeros(npar,1);
lbcum = zeros(npar,1);
ubcum = ones(npar,1);
bounds = prior_bounds(bayestopt_,options_);
[junk1,junk2,junk3,lb,ub,junk4] = set_prior(estim_params_,M_,options_); %Prepare bounds
if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
% Set prior bounds
bounds = prior_bounds(bayestopt_,options_);
bounds.lb = max(bounds.lb,lb);
bounds.ub = min(bounds.ub,ub);
else % estimated parameters but no declared priors
% No priors are declared so Dynare will estimate the model by
% maximum likelihood with inequality constraints for the parameters.
bounds.lb = lb;
bounds.ub = ub;
end
% set bounds for cumulative probabilities
for i = 1:npar
switch pshape(i)
......@@ -63,8 +74,8 @@ if init
p4(i) = min(p4(i),bounds.ub(i));
p3(i) = max(p3(i),bounds.lb(i));
case 3% Gaussian prior.
lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));;
ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));;
lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));
ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));
case 2% Gamma prior.
lbcum(i) = gamcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = gamcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
......
function set_shocks_param(xparam1)
% function set_shocks_param(xparam1)
% Set the structural and measurement error variances and covariances
% Copyright (C) 2012 Dynare Team
% Copyright (C) 2012-2015 Dynare Team
%
% This file is part of Dynare.
%
......@@ -18,31 +20,86 @@ function set_shocks_param(xparam1)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global estim_params_ M_
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
np = estim_params_.np;
Sigma_e = M_.Sigma_e;
offset = 0;
if nvx
offset = offset + nvx;
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
Sigma_e = M_.Sigma_e;
Correlation_matrix = M_.Correlation_matrix;
H = M_.H;
Correlation_matrix_ME = M_.Correlation_matrix_ME;
% setting shocks variance on the diagonal of Covariance matrix; used later
% for updating covariances
if nvx
var_exo = estim_params_.var_exo;
for i=1:nvx
k = var_exo(i,1);
Sigma_e(k,k) = xparam1(i)^2;
k =var_exo(i,1);
Sigma_e(k,k) = xparam1(i)^2;
end
end
if ncx
offset = offset + estim_params_.nvn;
end
% update offset
offset = nvx;
% setting measument error variance; on the diagonal of Covariance matrix; used later
% for updating covariances
if nvn
for i=1:nvn
k = estim_params_.nvn_observable_correspondence(i,1);
H(k,k) = xparam1(i+offset)^2;
end
end
% update offset
offset = nvx+nvn;
% setting shocks covariances
if ncx
corrx = estim_params_.corrx;
for i=1:ncx
k1 = corrx(i,1);
k2 = corrx(i,2);
Sigma_e(k1,k2) = xparam1(i+offset)*sqrt(Sigma_e_(k1,k1)*Sigma_e_(k2,k2));
Sigma_e(k2,k1) = Sigma_e_(k1,k2);
k1 = corrx(i,1);
k2 = corrx(i,2);
Correlation_matrix(k1,k2) = xparam1(i+offset);
Correlation_matrix(k2,k1) = Correlation_matrix(k1,k2);
end
end
end
%build covariance matrix from correlation matrix and variances already on
%diagonal
Sigma_e = diag(sqrt(diag(Sigma_e)))*Correlation_matrix*diag(sqrt(diag(Sigma_e)));
%if calibrated covariances, set them now to their stored value
if isfield(estim_params_,'calibrated_covariances')
Sigma_e(estim_params_.calibrated_covariances.position)=estim_params_.calibrated_covariances.cov_value;
end
% update offset
offset = nvx+nvn+ncx;
% setting measurement error covariances
if ncn
corrn_observable_correspondence = estim_params_.corrn_observable_correspondence;
for i=1:ncn
k1 = corrn_observable_correspondence(i,1);
k2 = corrn_observable_correspondence(i,2);
Correlation_matrix_ME(k1,k2) = xparam1(i+offset);
Correlation_matrix_ME(k2,k1) = Correlation_matrix_ME(k1,k2);
end
end
%build covariance matrix from correlation matrix and variances already on
%diagonal
H = diag(sqrt(diag(H)))*Correlation_matrix_ME*diag(sqrt(diag(H)));
%if calibrated covariances, set them now to their stored value
if isfield(estim_params_,'calibrated_covariances_ME')
H(estim_params_.calibrated_covariances_ME.position)=estim_params_.calibrated_covariances_ME.cov_value;
end
M_.Sigma_e = Sigma_e;
\ No newline at end of file
% updating matrices in M
if nvx || ncx
M_.Sigma_e = Sigma_e;
M_.Correlation_matrix=Correlation_matrix;
end
if nvn || ncn
M_.H = H;
M_.Correlation_matrix_ME=Correlation_matrix_ME;
end
\ No newline at end of file
......@@ -87,7 +87,7 @@ nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
lpmat0=[];
lpmat0=zeros(Nsam,0);
xparam1=[];
pshape = bayestopt_.pshape(nshock+1:end);
......@@ -96,7 +96,22 @@ p2 = bayestopt_.p2(nshock+1:end);
p3 = bayestopt_.p3(nshock+1:end);
p4 = bayestopt_.p4(nshock+1:end);
bounds = prior_bounds(bayestopt_,options_);
[junk1,junk2,junk3,lb,ub,junk4] = set_prior(estim_params_,M_,options_); %Prepare bounds
if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
% Set prior bounds
bounds = prior_bounds(bayestopt_,options_);
bounds.lb = max(bounds.lb,lb);
bounds.ub = min(bounds.ub,ub);
else % estimated parameters but no declared priors
% No priors are declared so Dynare will estimate the model by
% maximum likelihood with inequality constraints for the parameters.
bounds.lb = lb;
bounds.ub = ub;
if opt_gsa.prior_range==0
warning('GSA:: When using ML, sampling from the prior is not possible. Setting prior_range=1')
opt_gsa.prior_range=1;
end
end
if nargin==0,
OutputDirectoryName='';
......@@ -143,7 +158,7 @@ if fload==0,
end
end
else %ilptau==0
%[lpmat] = rand(Nsam,np);
[lpmat] = NaN(Nsam,np);
for j=1:np,
lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
end
......@@ -151,7 +166,7 @@ if fload==0,
end
end
% try
dummy=prior_draw_gsa(1);
dummy=prior_draw_gsa(1); %initialize persistent variables
% catch
% if pprior,
% if opt_gsa.prior_range==0;
......@@ -242,6 +257,7 @@ if fload==0,
else
d = chol(inv(hh));
lp=randn(Nsam*2,nshock+np)*d+kron(ones(Nsam*2,1),xparam1');
lnprior=zeros(1,Nsam*2);
for j=1:Nsam*2,
lnprior(j) = any(lp(j,:)'<=bounds.lb | lp(j,:)'>=bounds.ub);
end
......@@ -263,6 +279,7 @@ if fload==0,
iwrong=zeros(1,Nsam);
inorestriction=zeros(1,Nsam);
irestriction=zeros(1,Nsam);
infox=zeros(1,Nsam);
for j=1:Nsam,
M_ = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',estim_params_,M_);
%try stoch_simul([]);
......@@ -273,7 +290,7 @@ if fload==0,
[Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
end
infox(j,1)=info(1);
if infox(j,1)==0 && ~exist('T'),
if infox(j,1)==0 && ~exist('T','var'),
dr_=oo_.dr;
if prepSA,
try
......@@ -321,7 +338,7 @@ if fload==0,
% bayestopt_.restrict_columns, ...
% bayestopt_.restrict_aux);
end
if ~exist('nspred'),
if ~exist('nspred','var'),
nspred = dr_.nspred; %size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
......@@ -339,7 +356,7 @@ if fload==0,
istable(j)=0;
if isfield(dr_,'eigval')
egg(:,j) = sort(dr_.eigval);
if exist('nspred')
if exist('nspred','var')
if any(isnan(egg(1:nspred,j)))
iwrong(j)=j;
else
......@@ -349,7 +366,7 @@ if fload==0,
end
end
else
if exist('egg'),
if exist('egg','var'),
egg(:,j)=ones(size(egg,1),1).*NaN;
end
iwrong(j)=j;
......@@ -457,6 +474,7 @@ else
stoch_simul([]);
%T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable));
ntrans=length(istable);
yys=NaN(length(ys_),ntrans);
for j=1:ntrans,
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
%stoch_simul([]);
......@@ -465,12 +483,12 @@ else
%[Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
% bayestopt_.restrict_columns,...
% bayestopt_.restrict_aux);
if ~exist('T')
if ~exist('T','var')
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),ntrans);
end
dr_ = oo_.dr;
T(:,:,j) = [dr_.ghx dr_.ghu];
if ~exist('nspred')
if ~exist('nspred','var')
nspred = dr_.nspred; %size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
......
function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
% [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
% Copyright (C) 2012 Dynare Team
% Copyright (C) 2012-2015 Dynare Team
%
% This file is part of Dynare.
%
......@@ -17,38 +18,38 @@ function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ oo_ options_
nvar = size(var_list,1);
if nvar == 0
global M_ oo_ options_
nvar = size(var_list,2);
if nvar == 0
nvar = length(dr.order_var);
ivar = [1:nvar]';
else
else
ivar=zeros(nvar,1);
for i=1:nvar
i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
if isempty(i_tmp)
error (['One of the variable specified does not exist']) ;
else
ivar(i) = i_tmp;
end
i_tmp = strmatch(var_list(:,i),M_.endo_names,'exact');
if isempty(i_tmp)
error(['One of the variables specified does not exist']) ;
else
ivar(i) = i_tmp;
end
end
end
[gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_, options_);
m = dr.ys(ivar(stationary_vars));
end
[gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_, options_);
m = dr.ys(ivar(stationary_vars));
% i1 = find(abs(diag(gamma_y{1})) > 1e-12);
i1 = [1:length(ivar)];
s2 = diag(gamma_y{1});
sd = sqrt(s2);
i1 = [1:length(ivar)];
s2 = diag(gamma_y{1});
sd = sqrt(s2);
z = [ m sd s2 ];
mean = m;
var = gamma_y{1};
z = [ m sd s2 ];
mean = m;
var = gamma_y{1};
%'THEORETICAL MOMENTS';
%'MEAN','STD. DEV.','VARIANCE');
......@@ -56,28 +57,27 @@ z;
%'VARIANCE DECOMPOSITION (in percent)';
if M_.exo_nbr>1,
vdec = 100*gamma_y{options_.ar+2}(i1,:);
vdec = 100*gamma_y{options_.ar+2}(i1,:);
else
vdec = 100*ones(size(gamma_y{1}(i1,1)));
end
vdec = 100*ones(size(gamma_y{1}(i1,1)));
end
%'MATRIX OF CORRELATIONS';
if options_.opt_gsa.useautocorr,
corr = gamma_y{1}(i1,i1)./(sd(i1)*sd(i1)');
corr = corr-diag(diag(corr))+diag(diag(gamma_y{1}(i1,i1)));
else
corr = gamma_y{1}(i1,i1);
corr = gamma_y{1}(i1,i1);
end
if options_.ar > 0
%'COEFFICIENTS OF AUTOCORRELATION';
if options_.ar > 0
%'COEFFICIENTS OF AUTOCORRELATION';
for i=1:options_.ar
if options_.opt_gsa.useautocorr,
autocorr{i} = gamma_y{i+1}(i1,i1);
else
autocorr{i} = gamma_y{i+1}(i1,i1).*(sd(i1)*sd(i1)');
end
zz(:,i) = diag(gamma_y{i+1}(i1,i1));
if options_.opt_gsa.useautocorr,
autocorr{i} = gamma_y{i+1}(i1,i1);
else
autocorr{i} = gamma_y{i+1}(i1,i1).*(sd(i1)*sd(i1)');
end
zz(:,i) = diag(gamma_y{i+1}(i1,i1));
end
end
end
\ No newline at end of file
......@@ -86,7 +86,7 @@ if info(1)==0,
if init,
indJJ = (find(max(abs(JJ'),[],1)>1.e-8));
if isempty(indJJ) && any(any(isnan(JJ)))
error('There are NaN in the JJ matrix. Please check whether your model has units roots and you forgot to set lik_init~=1.' )
error('There are NaN in the JJ matrix. Please check whether your model has units roots and you forgot to set diffuse_filter=1.' )
end
while length(indJJ)<nparam && nlags<10,
disp('The number of moments with non-zero derivative is smaller than the number of parameters')
......@@ -118,14 +118,14 @@ if info(1)==0,
ide_strength_J_prior=NaN(1,nparam);
if init, %~isempty(indok),
normaliz = abs(params);
if prior_exist,
if prior_exist,
if ~isempty(estim_params_.var_exo),
normaliz1 = estim_params_.var_exo(:,7); % normalize with prior standard deviation
normaliz1 = estim_params_.var_exo(:,7)'; % normalize with prior standard deviation
else
normaliz1=[];
end
if ~isempty(estim_params_.param_vals),
normaliz1 = [normaliz1; estim_params_.param_vals(:,7)]'; % normalize with prior standard deviation
normaliz1 = [normaliz1 estim_params_.param_vals(:,7)']; % normalize with prior standard deviation
end
% normaliz = max([normaliz; normaliz1]);
normaliz1(isinf(normaliz1)) = 1;
......
......@@ -82,7 +82,7 @@ end
ixnoJ = 0;
if rankJ<npar || rankJJ<npar || min(1-McoJ)<1.e-10
if npar>0 && (rankJ<npar || rankJJ<npar || min(1-McoJ)<1.e-10)
% - find out which parameters are involved
% disp('Some parameters are NOT identified by the moments included in J')
% disp(' ')
......
......@@ -56,14 +56,22 @@ if SampleSize == 1,
subplot(211)
mmm = (idehess.ide_strength_J);
[ss, is] = sort(mmm);
bar(log([idehess.ide_strength_J(:,is)' idehess.ide_strength_J_prior(:,is)']))
if ~all(isnan(idehess.ide_strength_J_prior))
bar(log([idehess.ide_strength_J(:,is)' idehess.ide_strength_J_prior(:,is)']))
else
bar(log([idehess.ide_strength_J(:,is)' ]))
end
set(gca,'xlim',[0 nparam+1])
set(gca,'xticklabel','')
dy = get(gca,'ylim');
for ip=1:nparam,
text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
legend('relative to param value','relative to prior std','Location','Best')
if ~all(isnan(idehess.ide_strength_J_prior))