v4: sensitivity code reshaped and optimizer updated [Marco]
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@679 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
michel
parent
a8e01d23bf
commit
71cf74483f
|
|
@ -55,11 +55,6 @@ ifil2 = 1;
|
|||
ifil3 = 1;
|
||||
ifil4 = 1;
|
||||
h = waitbar(0,'Bayesian smoother...');
|
||||
if B >= MAX_nirfs
|
||||
stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,MAX_nirfs);
|
||||
else
|
||||
stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,B);
|
||||
end
|
||||
if options_.smoother
|
||||
if B <= MAX_nsmoo
|
||||
stock_smooth = zeros(endo_nbr,gend,B);
|
||||
|
|
@ -106,38 +101,42 @@ for b=1:B
|
|||
irun3 = irun3 + 1;
|
||||
irun4 = irun4 + 1;
|
||||
|
||||
if irun1 == MAX_nsmoo | b == B
|
||||
if irun1 > MAX_nsmoo | b == B
|
||||
if b == B
|
||||
stock_smooth = stock_smoo(:,:,1:irun1);
|
||||
end
|
||||
save([DirectoryName M_.fname '_smooth' int2str(ifil1)],'stock_smooth');
|
||||
stock = stock_smooth;
|
||||
save([DirectoryName M_.fname '_smooth' int2str(ifil1)],'stock');
|
||||
ifil1 = ifil1 + 1;
|
||||
irun1 = 1;
|
||||
end
|
||||
|
||||
if nvx & (irun2 == MAX_inno | b == B)
|
||||
if nvx & (irun2 > MAX_inno | b == B)
|
||||
if b == B
|
||||
stock_innov = stock_innov(:,:,1:irun2);
|
||||
end
|
||||
save([DirectoryName M_.fname '_inno' int2str(ifil2)],'stock_inno');
|
||||
stock = stock_inno'
|
||||
save([DirectoryName M_.fname '_inno' int2str(ifil2)],'stock');
|
||||
ifil2 = ifil2 + 1;
|
||||
irun2 = 1;
|
||||
end
|
||||
|
||||
if nvn & (irun3 == MAX_error | b == B)
|
||||
if nvn & (irun3 > MAX_error | b == B)
|
||||
if b == B
|
||||
stock_error = stock_error(:,:,1:irun3);
|
||||
end
|
||||
save([DirectoryName M_.fname '_error' int2str(ifil3)],'stock_error');
|
||||
stock = stock_error;
|
||||
save([DirectoryName M_.fname '_error' int2str(ifil3)],'stock');
|
||||
ifil3 = ifil3 + 1;
|
||||
irun3 = 1;
|
||||
end
|
||||
|
||||
if naK & (irun3 == MAX_naK | b == B)
|
||||
if naK & (irun3 > MAX_naK | b == B)
|
||||
if b == B
|
||||
stock_filter = stock_filter(:,:,:,1:irun4);
|
||||
end
|
||||
save([DirectoryName M_.fname '_filter' int2str(ifil4)],'stock_filter');
|
||||
stock = stock_filter;
|
||||
save([DirectoryName M_.fname '_filter' int2str(ifil4)],'stock');
|
||||
ifil4 = ifil4 + 1;
|
||||
irun4 = 1;
|
||||
end
|
||||
|
|
@ -145,3 +144,4 @@ for b=1:B
|
|||
waitbar(b/B,h);
|
||||
end
|
||||
close(h)
|
||||
|
||||
|
|
|
|||
|
|
@ -3,6 +3,9 @@ function dynare_estimation(var_list_)
|
|||
global M_ options_ oo_ estim_params_
|
||||
global bayestopt_
|
||||
|
||||
% temporary fix until M_.H is initialized by the parser
|
||||
M_.H = [];
|
||||
|
||||
options_.varlist = var_list_;
|
||||
options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1);
|
||||
for i = 1:size(M_.endo_names,1)
|
||||
|
|
@ -58,6 +61,7 @@ options_ = set_default_option(options_,'MaxNumberOfBytes',1e6);
|
|||
options_ = set_default_option(options_,'xls_sheet','');
|
||||
options_ = set_default_option(options_,'xls_range','');
|
||||
options_ = set_default_option(options_,'filter_step_ahead',0);
|
||||
options_ = set_default_option(options_,'diffuse_d',[]);
|
||||
if options_.filtered_vars ~= 0 & options_filter_step_ahead == 0
|
||||
options_.filter_step_ahead = 1;
|
||||
end
|
||||
|
|
@ -256,9 +260,29 @@ if options_.mode_compute > 0 & options_.posterior_mode_estimation
|
|||
disp(sprintf('Objective function at mode: %f',fval))
|
||||
disp(sprintf('Objective function at mode: %f',DsgeLikelihood(xparam1,gend,data)))
|
||||
elseif options_.mode_compute == 5
|
||||
flag = 0;
|
||||
[xparam1, hh, gg, fval] = newrat('DsgeLikelihood',xparam1,[],[],flag,gend,data);
|
||||
eval(['save ' M_.fname '_mode xparam1 hh gg fval;']);
|
||||
if isfield(options_,'hess')
|
||||
flag = options_.hess;
|
||||
else
|
||||
flag = 1;
|
||||
end
|
||||
if ~exist('igg'), % by M. Ratto
|
||||
hh=[];
|
||||
gg=[];
|
||||
igg=[];
|
||||
end % by M. Ratto
|
||||
if isfield(options_,'ftol')
|
||||
crit = options_.ftol;
|
||||
else
|
||||
crit = 1.e-7;
|
||||
end
|
||||
if isfield(options_,'nit')
|
||||
nit = options_.nit;
|
||||
else
|
||||
nit=1000;
|
||||
end
|
||||
%[xparam1, hh, gg, fval] = newrat('DsgeLikelihood',xparam1,[],[],flag,gend,data);
|
||||
[xparam1, hh, gg, fval, invhess] = newrat('DsgeLikelihood',xparam1,hh,gg,igg,crit,nit,flag,gend,data);
|
||||
eval(['save ' M_.fname '_mode xparam1 hh gg fval invhess;']);
|
||||
end
|
||||
if options_.mode_compute ~= 5
|
||||
hh = reshape(hessian('DsgeLikelihood',xparam1,gend,data),nx,nx);
|
||||
|
|
@ -665,9 +689,8 @@ if (any(bayestopt_.pshape >0 ) & options_.mh_replic) | ...
|
|||
%%
|
||||
GetPosteriorParametersStatistics;
|
||||
PlotPosteriorDistributions;
|
||||
metropolis_draw(1);% Initialization.
|
||||
if options_.bayesian_irf
|
||||
PosteriorIRF('posterior');
|
||||
PosteriorIRF;
|
||||
end
|
||||
return
|
||||
|
||||
|
|
|
|||
|
|
@ -34,11 +34,13 @@ icount=0;
|
|||
nx=length(x);
|
||||
xparam1=x;
|
||||
%ftol0=1.e-6;
|
||||
htol_base = max(1.e-5, ftol0);
|
||||
flagit=0; % mode of computation of hessian in each iteration
|
||||
ftol=ftol0;
|
||||
gtol=1.e-3;
|
||||
htol=ftol0;
|
||||
htol0=ftol0;
|
||||
htol=htol_base;
|
||||
htol0=htol_base;
|
||||
gibbstol=length(bayestopt_.pshape)/12;
|
||||
|
||||
func_hh = [func0,'_hh'];
|
||||
func = str2func(func0);
|
||||
|
|
@ -56,7 +58,7 @@ if isempty(hh)
|
|||
end
|
||||
if htol0>htol,
|
||||
htol=htol0;
|
||||
ftol=htol0;
|
||||
%ftol=htol0;
|
||||
end
|
||||
else
|
||||
hh0=hh;
|
||||
|
|
@ -85,18 +87,18 @@ while norm(gg)>gtol & check==0 & jit<nit,
|
|||
disp(['Iteration ',num2str(icount)])
|
||||
[fval x0 fc retcode] = csminit(func0,xparam1,fval0(icount),gg,0,igg,varargin{:});
|
||||
if igrad,
|
||||
[fval1 x01 fc retcode1] = csminit(func0,xparam1,fval0(icount),gg,0,inx,varargin{:});
|
||||
if fval1<fval,
|
||||
fval=fval1;
|
||||
x0=x01;
|
||||
[fval1 x01 fc retcode1] = csminit(func0,x0,fval,gg,0,inx,varargin{:});
|
||||
if (fval-fval1)>1, %(fval0(icount)-fval),
|
||||
disp('Gradient step!!')
|
||||
else
|
||||
igrad=0;
|
||||
end
|
||||
fval=fval1;
|
||||
x0=x01;
|
||||
end
|
||||
if (fval0(icount)-fval)<1.e-2*(gg'*(igg*gg))/2 & igibbs,
|
||||
[fvala, x0] = mr_gstep(func0,x0,htol,varargin{:});
|
||||
if (fval-fvala)<5*(fval0(icount)-fval),
|
||||
if (fval-fvala)<gibbstol*(fval0(icount)-fval),
|
||||
igibbs=0;
|
||||
disp('Last Gibbs step, gain too small!!')
|
||||
else
|
||||
|
|
@ -107,29 +109,23 @@ while norm(gg)>gtol & check==0 & jit<nit,
|
|||
if (fval0(icount)-fval)<ftol & flagit==0,
|
||||
disp('Try diagonal Hessian')
|
||||
ihh=diag(1./(diag(hhg)));
|
||||
[fval2 x02 fc retcode2] = csminit(func2str(func),xparam1,fval0(icount),gg,0,ihh,varargin{:});
|
||||
if fval2<fval,
|
||||
x0=x02;
|
||||
fval=fval2;
|
||||
if (fval0(icount)-fval2)>=ftol ,
|
||||
[fval2 x0 fc retcode2] = csminit(func2str(func),x0,fval,gg,0,ihh,varargin{:});
|
||||
if (fval-fval2)>=ftol ,
|
||||
%hh=diag(diag(hh));
|
||||
disp('Diagonal Hessian successful')
|
||||
end
|
||||
end
|
||||
fval=fval2;
|
||||
end
|
||||
if (fval0(icount)-fval)<ftol & flagit==0,
|
||||
disp('Try gradient direction')
|
||||
ihh0=inx.*1.e-4;
|
||||
[fval3 x03 fc retcode3] = csminit(func2str(func),xparam1,fval0(icount),gg,0,ihh0,varargin{:});
|
||||
if fval3<fval,
|
||||
x0=x03;
|
||||
fval=fval3;
|
||||
if (fval0(icount)-fval3)>=ftol ,
|
||||
[fval3 x0 fc retcode3] = csminit(func2str(func),x0,fval,gg,0,ihh0,varargin{:});
|
||||
if (fval-fval3)>=ftol ,
|
||||
%hh=hh0;
|
||||
%ihh=ihh0;
|
||||
disp('Gradient direction successful')
|
||||
end
|
||||
end
|
||||
fval=fval3;
|
||||
end
|
||||
xparam1=x0;
|
||||
x(:,icount+1)=xparam1;
|
||||
|
|
@ -172,21 +168,21 @@ while norm(gg)>gtol & check==0 & jit<nit,
|
|||
disp(['Htol ',num2str(htol0)])
|
||||
|
||||
if df<htol0,
|
||||
htol=max(ftol0,df/10);
|
||||
htol=max(htol_base,df/10);
|
||||
end
|
||||
|
||||
if norm(x(:,icount)-xparam1)>1.e-12,
|
||||
save m1 x fval0 -append
|
||||
[dum, gg, htol0, igg, hhg]=mr_hessian(func_hh,xparam1,flagit,htol,varargin{:});
|
||||
if htol0>ftol,
|
||||
ftol=htol0;
|
||||
if htol0>htol, %ftol,
|
||||
%ftol=htol0;
|
||||
htol=htol0;
|
||||
disp(' ')
|
||||
disp('Numerical noise in the likelihood')
|
||||
disp('Tolerance has to be relaxed')
|
||||
disp(' ')
|
||||
elseif htol0<ftol,
|
||||
ftol=max(htol0, ftol0);
|
||||
% elseif htol0<ftol,
|
||||
% ftol=max(htol0, ftol0);
|
||||
end
|
||||
hh0 = reshape(dum,nx,nx);
|
||||
hh=hhg;
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
function [H,prob,d] = smirnov(x1 , x2 , alpha )
|
||||
function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
|
||||
% Smirnov test for 2 distributions
|
||||
% [H,prob,d] = smirnov(x1 , x2 , alpha )
|
||||
% [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
|
||||
%
|
||||
% Copyright (C) 2005 Marco Ratto
|
||||
%
|
||||
|
|
@ -10,40 +10,49 @@ function [H,prob,d] = smirnov(x1 , x2 , alpha )
|
|||
if nargin<3
|
||||
alpha = 0.05;
|
||||
end
|
||||
|
||||
if nargin<4,
|
||||
iflag=0;
|
||||
end
|
||||
|
||||
% empirical cdfs.
|
||||
xmix= [x1;x2];
|
||||
bin = [-inf ; sort(xmix) ; inf];
|
||||
|
||||
bin = [-inf ; sort([x1;x2]) ; inf];
|
||||
|
||||
ncount1 = histc (x1 , bin);
|
||||
ncount2 = histc (x2 , bin);
|
||||
ncount1 = histc (x1 , bin);
|
||||
ncount2 = histc (x2 , bin);
|
||||
|
||||
cum1 = cumsum(ncount1)./sum(ncount1);
|
||||
cum2 = cumsum(ncount2)./sum(ncount2);
|
||||
|
||||
cum1 = cum1(1:end-1);
|
||||
|
||||
cum2 = cumsum(ncount2)./sum(ncount2);
|
||||
cum2 = cum2(1:end-1);
|
||||
|
||||
n1 = length(x1);
|
||||
n2 = length(x2);
|
||||
n = n1 * n2 /(n1 + n2);
|
||||
n1= length(x1);
|
||||
n2= length(x2);
|
||||
n = n1*n2 /(n1+n2);
|
||||
|
||||
% Compute the d(n1,n2) statistic.
|
||||
|
||||
d = max(abs(cum1 - cum2));
|
||||
% Compute the d(n1,n2) statistics.
|
||||
|
||||
if iflag==0,
|
||||
d = max(abs(cum1 - cum2));
|
||||
elseif iflag==-1
|
||||
d = max(cum2 - cum1);
|
||||
elseif iflag==1
|
||||
d = max(cum1 - cum2);
|
||||
end
|
||||
%
|
||||
% Compute P-value check H0 hypothesis
|
||||
%
|
||||
|
||||
lam = max((sqrt(n) + 0.12 + 0.11/sqrt(n)) * d , 0);
|
||||
|
||||
j = [1:101]';
|
||||
prob = 2 * sum((-1).^(j-1).*exp(-2*lam*lam*j.^2));
|
||||
|
||||
if prob < 0 , prob = 0; end
|
||||
if prob > 1 , prob = 1; end
|
||||
|
||||
if iflag == 0
|
||||
j = [1:101]';
|
||||
prob = 2 * sum((-1).^(j-1).*exp(-2*lam*lam*j.^2));
|
||||
|
||||
prob=max(prob,0);
|
||||
prob=min(1,prob);
|
||||
else
|
||||
prob = exp(-2*lam*lam);
|
||||
end
|
||||
|
||||
H = (alpha >= prob);
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
function x0 = stab_map_(Nsam, fload, alpha2, alpha)
|
||||
function x0 = stab_map_(Nsam, fload, alpha2, alpha, prepSA, pprior)
|
||||
%
|
||||
% function x0 = stab_map_(Nsam, fload, alpha2, alpha)
|
||||
% function x0 = stab_map_(Nsam, fload, alpha2, alpha, prepSA, pprior)
|
||||
%
|
||||
% Mapping of stability regions in the prior ranges applying
|
||||
% Monte Carlo filtering techniques.
|
||||
|
|
@ -15,7 +15,12 @@ function x0 = stab_map_(Nsam, fload, alpha2, alpha)
|
|||
% [abs(corrcoef) > alpha2]
|
||||
% alpha = significance level for univariate sensitivity analysis
|
||||
% (uses smirnov)
|
||||
%
|
||||
% prepSA = 1: save transition matrices for mapping reduced form
|
||||
% = 0: no transition matrix saved (default)
|
||||
% pprior = 1: sample from prior ranges (default): sample saved in
|
||||
% _stab.mat file
|
||||
% = 0: sample from posterior ranges: sample saved in
|
||||
% _mc.mat file
|
||||
% OUTPUT:
|
||||
% x0: one parameter vector for which the model is stable.
|
||||
%
|
||||
|
|
@ -29,7 +34,7 @@ function x0 = stab_map_(Nsam, fload, alpha2, alpha)
|
|||
% ( abs(corrcoef) > alpha2) under the stable subset
|
||||
%
|
||||
% USES lptauSEQ,
|
||||
% smirnov
|
||||
% stab_map_1, stab_map_2
|
||||
%
|
||||
% Copyright (C) 2005 Marco Ratto
|
||||
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
|
||||
|
|
@ -61,10 +66,6 @@ nshock = estim_params_.nvx;
|
|||
nshock = nshock + estim_params_.nvn;
|
||||
nshock = nshock + estim_params_.ncx;
|
||||
nshock = nshock + estim_params_.ncn;
|
||||
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
|
||||
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
|
||||
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
|
||||
%kernel_function = 'uniform'; % Gaussian kernel for Fast Fourrier Transform approximaton.
|
||||
|
||||
|
||||
if nargin==0,
|
||||
|
|
@ -73,13 +74,24 @@ end
|
|||
if nargin<2,
|
||||
fload=0;
|
||||
end
|
||||
if nargin<4,
|
||||
alpha=0.002;
|
||||
end
|
||||
if nargin<3,
|
||||
alpha2=0.3;
|
||||
end
|
||||
if nargin<4 | isempty(alpha),
|
||||
alpha=0.002;
|
||||
end
|
||||
if nargin<5,
|
||||
prepSA=0;
|
||||
end
|
||||
if nargin<6,
|
||||
pprior=1;
|
||||
end
|
||||
|
||||
options_.periods=0;
|
||||
options_.nomoments=1;
|
||||
options_.irf=0;
|
||||
options_.noprint=1;
|
||||
|
||||
if fload==0 | nargin<2 | isempty(fload),
|
||||
if estim_params_.np<52,
|
||||
[lpmat] = lptauSEQ(Nsam,estim_params_.np);
|
||||
|
|
@ -90,39 +102,71 @@ if fload==0 | nargin<2 | isempty(fload),
|
|||
end
|
||||
end
|
||||
|
||||
|
||||
for j=1:estim_params_.np,
|
||||
if estim_params_.np>30 & estim_params_.np<52
|
||||
lpmat(:,j)=lpmat(randperm(Nsam),j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
|
||||
else
|
||||
lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
|
||||
if pprior,
|
||||
for j=1:estim_params_.np,
|
||||
if estim_params_.np>30 & estim_params_.np<52
|
||||
lpmat(:,j)=lpmat(randperm(Nsam),j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
|
||||
else
|
||||
lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
|
||||
end
|
||||
end
|
||||
else
|
||||
for j=1:nshock,
|
||||
xparam1(j) = oo_.posterior_mode.shocks_std.(bayestopt_.name{j});
|
||||
sd(j) = oo_.posterior_std.shocks_std.(bayestopt_.name{j});
|
||||
lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
|
||||
lb = max(bayestopt_.lb(j), xparam1(j)-2*sd(j));
|
||||
ub1=xparam1(j)+(xparam1(j) - lb); % define symmetric range around the mode!
|
||||
ub = min(bayestopt_.ub(j),ub1);
|
||||
if ub<ub1,
|
||||
lb=xparam1(j)-(ub-xparam1(j)); % define symmetric range around the mode!
|
||||
end
|
||||
lpmat0(:,j) = lpmat0(:,j).*(ub-lb)+lb;
|
||||
end
|
||||
%
|
||||
for j=1:estim_params_.np,
|
||||
xparam1(j+nshock) = oo_.posterior_mode.parameters.(bayestopt_.name{j+nshock});
|
||||
sd(j+nshock) = oo_.posterior_std.parameters.(bayestopt_.name{j+nshock});
|
||||
lb = max(bayestopt_.lb(j+nshock),xparam1(j+nshock)-2*sd(j+nshock));
|
||||
ub1=xparam1(j+nshock)+(xparam1(j+nshock) - lb); % define symmetric range around the mode!
|
||||
ub = min(bayestopt_.ub(j+nshock),ub1);
|
||||
if ub<ub1,
|
||||
lb=xparam1(j+nshock)-(ub-xparam1(j+nshock)); % define symmetric range around the mode!
|
||||
end
|
||||
%ub = min(bayestopt_.ub(j+nshock),xparam1(j+nshock)+2*sd(j+nshock));
|
||||
if estim_params_.np>30 & estim_params_.np<52
|
||||
lpmat(:,j) = lpmat(randperm(Nsam),j).*(ub-lb)+lb;
|
||||
else
|
||||
lpmat(:,j) = lpmat(:,j).*(ub-lb)+lb;
|
||||
end
|
||||
end
|
||||
end
|
||||
%
|
||||
h = waitbar(0,'Please wait...');
|
||||
options_.periods=0;
|
||||
options_.nomoments=1;
|
||||
options_.irf=0;
|
||||
options_.noprint=1;
|
||||
istable=[1:Nsam];
|
||||
iunstable=[1:Nsam];
|
||||
for j=1:Nsam,
|
||||
% for i=1:estim_params_.np,
|
||||
% evalin('base',[bayestopt_.name{i+nshock}, '= ',sprintf('%0.15e',lpmat(j,i)),';'])
|
||||
% end
|
||||
M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)';
|
||||
%evalin('base','stoch_simul(var_list_);');
|
||||
stoch_simul([]);
|
||||
dr_ = oo_.dr;
|
||||
%egg(:,j) = sort(eigenvalues_);
|
||||
%egg(:,j) = sort(dr_.eigval);
|
||||
if isfield(dr_,'eigval'),
|
||||
if isfield(dr_,'ghx'),
|
||||
egg(:,j) = sort(dr_.eigval);
|
||||
if ~exist('nspred')
|
||||
iunstable(j)=0;
|
||||
if prepSA
|
||||
T(:,:,j) = [dr_.ghx dr_.ghu];
|
||||
end
|
||||
if ~exist('nspred'),
|
||||
nspred = size(dr_.ghx,2);
|
||||
nboth = dr_.nboth;
|
||||
nfwrd = dr_.nfwrd;
|
||||
end
|
||||
else
|
||||
egg(:,j)=ones(size(egg,1),1).*1.1;
|
||||
istable(j)=0;
|
||||
if isfield(dr_,'eigval')
|
||||
egg(:,j) = sort(dr_.eigval);
|
||||
else
|
||||
egg(:,j)=ones(size(egg,1),1).*1.1;
|
||||
end
|
||||
end
|
||||
ys_=real(dr_.ys);
|
||||
yys(:,j) = ys_;
|
||||
|
|
@ -130,160 +174,140 @@ if fload==0 | nargin<2 | isempty(fload),
|
|||
waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
|
||||
end
|
||||
close(h)
|
||||
istable=istable(find(istable)); % stable params
|
||||
iunstable=iunstable(find(iunstable)); % unstable params
|
||||
|
||||
% map stable samples
|
||||
ix=[1:Nsam];
|
||||
for j=1:Nsam,
|
||||
if abs(egg(nspred,j))>=options_.qz_criterium; %(1-(options_.qz_criterium-1)); %1-1.e-5;
|
||||
ix(j)=0;
|
||||
%elseif (dr_.nboth | dr_.nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
|
||||
elseif (nboth | nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
|
||||
ix(j)=0;
|
||||
end
|
||||
end
|
||||
ix=ix(find(ix)); % stable params
|
||||
|
||||
% map unstable samples
|
||||
ixx=[1:Nsam];
|
||||
for j=1:Nsam,
|
||||
%if abs(egg(dr_.npred+1,j))>1+1.e-5 & abs(egg(dr_.npred,j))<1-1.e-5;
|
||||
%if (dr_.nboth | dr_.nfwrd),
|
||||
if (nboth | nfwrd),
|
||||
if abs(egg(nspred+1,j))>options_.qz_criterium & abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
|
||||
ixx(j)=0;
|
||||
end
|
||||
% % map stable samples
|
||||
% istable=[1:Nsam];
|
||||
% for j=1:Nsam,
|
||||
% if any(isnan(egg(1:nspred,j)))
|
||||
% istable(j)=0;
|
||||
% else
|
||||
% if abs(egg(nspred,j))>=options_.qz_criterium; %(1-(options_.qz_criterium-1)); %1-1.e-5;
|
||||
% istable(j)=0;
|
||||
% %elseif (dr_.nboth | dr_.nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
|
||||
% elseif (nboth | nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
|
||||
% istable(j)=0;
|
||||
% end
|
||||
% end
|
||||
% end
|
||||
% istable=istable(find(istable)); % stable params
|
||||
%
|
||||
% % map unstable samples
|
||||
% iunstable=[1:Nsam];
|
||||
% for j=1:Nsam,
|
||||
% %if abs(egg(dr_.npred+1,j))>1+1.e-5 & abs(egg(dr_.npred,j))<1-1.e-5;
|
||||
% %if (dr_.nboth | dr_.nfwrd),
|
||||
% if ~any(isnan(egg(1:5,j)))
|
||||
% if (nboth | nfwrd),
|
||||
% if abs(egg(nspred+1,j))>options_.qz_criterium & abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
|
||||
% iunstable(j)=0;
|
||||
% end
|
||||
% else
|
||||
% if abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
|
||||
% iunstable(j)=0;
|
||||
% end
|
||||
% end
|
||||
% end
|
||||
% end
|
||||
% iunstable=iunstable(find(iunstable)); % unstable params
|
||||
if pprior,
|
||||
if ~prepSA
|
||||
save([fname_,'_stab'],'lpmat','iunstable','istable','egg','yys','nspred','nboth','nfwrd')
|
||||
else
|
||||
if abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
|
||||
ixx(j)=0;
|
||||
end
|
||||
save([fname_,'_stab'],'lpmat','iunstable','istable','egg','yys','T','nspred','nboth','nfwrd')
|
||||
end
|
||||
end
|
||||
ixx=ixx(find(ixx)); % unstable params
|
||||
save([fname_,'_stab'],'lpmat','ixx','ix','egg','yys')
|
||||
|
||||
else
|
||||
if ~prepSA
|
||||
save([fname_,'_mc'],'lpmat','lpmat0','iunstable','istable','egg','yys','nspred','nboth','nfwrd')
|
||||
else
|
||||
save([fname_,'_mc'],'lpmat','lpmat0','iunstable','istable','egg','yys','T','nspred','nboth','nfwrd')
|
||||
end
|
||||
end
|
||||
else
|
||||
load([fname_,'_stab'])
|
||||
if pprior,
|
||||
load([fname_,'_stab'])
|
||||
else
|
||||
load([fname_,'_mc'])
|
||||
end
|
||||
Nsam = size(lpmat,1);
|
||||
end
|
||||
|
||||
delete([fname_,'_stab_*.*']);
|
||||
delete([fname_,'_stab_SA_*.*']);
|
||||
delete([fname_,'_stab_corr_*.*']);
|
||||
delete([fname_,'_unstab_corr_*.*']);
|
||||
|
||||
if length(ixx)>0 & length(ixx)<Nsam,
|
||||
% Blanchard Kahn
|
||||
for i=1:ceil(estim_params_.np/12),
|
||||
figure,
|
||||
for j=1+12*(i-1):min(estim_params_.np,12*i),
|
||||
subplot(3,4,j-12*(i-1))
|
||||
optimal_bandwidth = mh_optimal_bandwidth(lpmat(ix,j),length(ix),bandwidth,kernel_function);
|
||||
[x1,f1] = kernel_density_estimate(lpmat(ix,j),number_of_grid_points,...
|
||||
optimal_bandwidth,kernel_function);
|
||||
plot(x1, f1,':k','linewidth',2)
|
||||
optimal_bandwidth = mh_optimal_bandwidth(lpmat(ixx,j),length(ixx),bandwidth,kernel_function);
|
||||
[x1,f1] = kernel_density_estimate(lpmat(ixx,j),number_of_grid_points,...
|
||||
optimal_bandwidth,kernel_function);
|
||||
hold on, plot(x1, f1,'k','linewidth',2)
|
||||
|
||||
%hist(lpmat(ix,j),30)
|
||||
title(bayestopt_.name{j+nshock})
|
||||
end
|
||||
saveas(gcf,[fname_,'_stab_',int2str(i)])
|
||||
end
|
||||
|
||||
% Smirnov test for Blanchard;
|
||||
for i=1:ceil(estim_params_.np/12),
|
||||
figure,
|
||||
for j=1+12*(i-1):min(estim_params_.np,12*i),
|
||||
subplot(3,4,j-12*(i-1))
|
||||
if ~isempty(ix),
|
||||
h=cumplot(lpmat(ix,j));
|
||||
set(h,'color',[0 0 0], 'linestyle',':')
|
||||
if prepSA & ~exist('T'),
|
||||
h = waitbar(0,'Please wait...');
|
||||
options_.periods=0;
|
||||
options_.nomoments=1;
|
||||
options_.irf=0;
|
||||
options_.noprint=1;
|
||||
stoch_simul([]);
|
||||
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable));
|
||||
|
||||
for j=1:length(istable),
|
||||
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
|
||||
stoch_simul([]);
|
||||
dr_ = oo_.dr;
|
||||
T(:,:,j) = [dr_.ghx dr_.ghu];
|
||||
if ~exist('nspred')
|
||||
nspred = size(dr_.ghx,2);
|
||||
nboth = dr_.nboth;
|
||||
nfwrd = dr_.nfwrd;
|
||||
end
|
||||
hold on,
|
||||
if ~isempty(ixx),
|
||||
h=cumplot(lpmat(ixx,j));
|
||||
set(h,'color',[0 0 0])
|
||||
end
|
||||
% if exist('kstest2')==2 & length(ixx)>0 & length(ixx)<Nsam,
|
||||
% [H,P,KSSTAT] = kstest2(lpmat(ix,j),lpmat(ixx,j));
|
||||
% title([bayestopt_.name{j+nshock},'. K-S prob ', num2str(P)])
|
||||
% else
|
||||
[H,P,KSSTAT] = smirnov(lpmat(ix,j),lpmat(ixx,j));
|
||||
title([bayestopt_.name{j+nshock},'. K-S prob ', num2str(P)])
|
||||
% end
|
||||
ys_=real(dr_.ys);
|
||||
yys(:,j) = ys_;
|
||||
ys_=yys(:,1);
|
||||
waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
|
||||
end
|
||||
saveas(gcf,[fname_,'_stab_SA_',int2str(i)])
|
||||
end
|
||||
|
||||
|
||||
disp(' ')
|
||||
disp(' ')
|
||||
disp('Starting bivariate analysis:')
|
||||
|
||||
c0=corrcoef(lpmat(ix,:));
|
||||
c00=tril(c0,-1);
|
||||
|
||||
stab_map_2(lpmat(ix,:),alpha2, 1);
|
||||
stab_map_2(lpmat(ixx,:),alpha2, 0);
|
||||
|
||||
else
|
||||
if length(ixx)==0,
|
||||
disp('All parameter values in the prior ranges are stable!')
|
||||
close(h)
|
||||
if pprior
|
||||
save([fname_,'_stab'],'T','-append')
|
||||
else
|
||||
disp('All parameter values in the prior ranges are unstable!')
|
||||
save([fname_,'_mc'],'T','-append')
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
if pprior
|
||||
aname='stab';
|
||||
auname='unstable';
|
||||
asname='stable';
|
||||
else
|
||||
aname='mc_stab';
|
||||
auname='mc_unstable';
|
||||
asname='mc_stable';
|
||||
end
|
||||
delete([fname_,'_',aname,'_*.*']);
|
||||
delete([fname_,'_',aname,'_SA_*.*']);
|
||||
delete([fname_,'_',asname,'_corr_*.*']);
|
||||
delete([fname_,'_',auname,'_corr_*.*']);
|
||||
|
||||
% % optional map cyclicity of dominant eigenvalues, if
|
||||
% thex=[];
|
||||
% for j=1:Nsam,
|
||||
% %cyc(j)=max(abs(imag(egg(1:34,j))));
|
||||
% ic = find(imag(egg(1:nspred,j)));
|
||||
% i=find( abs(egg( ic ,j) )>0.9); %only consider complex dominant eigenvalues
|
||||
% if ~isempty(i),
|
||||
% i=i(1:2:end);
|
||||
% thedum=[];
|
||||
% for ii=1:length(i),
|
||||
% idum = ic( i(ii) );
|
||||
% thedum(ii)=abs(angle(egg(idum,j)));
|
||||
% end
|
||||
% [dum, icx]=max(thedum);
|
||||
% icy(j) = ic( i(icx) );
|
||||
% thet(j)=max(thedum);
|
||||
% if thet(j)<0.05 & find(ix==j), % keep stable runs with freq smaller than 0.05
|
||||
% thex=[thex; j];
|
||||
% end
|
||||
% else
|
||||
% if find(ix==j),
|
||||
% thex=[thex; j];
|
||||
% end
|
||||
% end
|
||||
% end
|
||||
% % cyclicity
|
||||
% for i=1:ceil(estim_params_.np/12),
|
||||
% figure,
|
||||
% for j=1+12*(i-1):min(estim_params_.np,12*i),
|
||||
% subplot(3,4,j-12*(i-1))
|
||||
% hist(lpmat(thex,j),30)
|
||||
% title(bayestopt_.name{j+nshock})
|
||||
% end
|
||||
% end
|
||||
%
|
||||
% % TFP STEP & Blanchard; & cyclicity
|
||||
% for i=1:ceil(estim_params_.np/12),
|
||||
% figure,
|
||||
% for j=1+12*(i-1):min(estim_params_.np,12*i),
|
||||
% [H,P,KSSTAT] = kstest2(lpmat(1:Nsam,j),lpmat(ixx,j));
|
||||
% subplot(3,4,j-12*(i-1))
|
||||
% cdfplot(lpmat(1:Nsam,j))
|
||||
% hold on,
|
||||
% cdfplot(lpmat(ixx,j))
|
||||
% title([bayestopt_.name{j+nshock},'. K-S prob ', num2str(P)])
|
||||
% end
|
||||
% end
|
||||
|
||||
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
|
||||
x0 = [x0; lpmat(ix(1),:)'];
|
||||
if length(iunstable)>0 & length(iunstable)<Nsam,
|
||||
disp([num2str(length(istable)/Nsam*100),'\% of the prior support is stable.'])
|
||||
% Blanchard Kahn
|
||||
proba = stab_map_1(lpmat, istable, iunstable, aname);
|
||||
disp(' ')
|
||||
disp(' ')
|
||||
disp('Starting bivariate analysis:')
|
||||
|
||||
c0=corrcoef(lpmat(istable,:));
|
||||
c00=tril(c0,-1);
|
||||
|
||||
stab_map_2(lpmat(istable,:),alpha2, asname);
|
||||
stab_map_2(lpmat(iunstable,:),alpha2, auname);
|
||||
|
||||
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
|
||||
x0 = [x0; lpmat(istable(1),:)'];
|
||||
if istable(end)~=Nsam
|
||||
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(1),:)';
|
||||
stoch_simul([]);
|
||||
end
|
||||
else
|
||||
if length(iunstable)==0,
|
||||
disp('All parameter values in the specified ranges are stable!')
|
||||
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
|
||||
x0 = [x0; lpmat(istable(1),:)'];
|
||||
else
|
||||
disp('All parameter values in the specified ranges are unstable!')
|
||||
x0=[];
|
||||
end
|
||||
|
||||
end
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
function stab_map_2(x,alpha2,istab)
|
||||
% function stab_map_2(x,alpha2,istab)
|
||||
%function stab_map_2(x,alpha2,istab,fnam)
|
||||
function stab_map_2(x,alpha2,fnam, ishock)
|
||||
% function stab_map_2(x,alpha2,fnam, ishock)
|
||||
%
|
||||
% Copyright (C) 2005 Marco Ratto
|
||||
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
|
||||
|
|
@ -12,19 +13,31 @@ function stab_map_2(x,alpha2,istab)
|
|||
% marco.ratto@jrc.it
|
||||
%
|
||||
|
||||
global bayestopt_ estim_params_ dr_ options_ ys_ fname_
|
||||
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
|
||||
global bayestopt_ estim_params_ options_ oo_ M_
|
||||
|
||||
if nargin<4,
|
||||
ishock=0;
|
||||
end
|
||||
ys_ = oo_.dr.ys;
|
||||
dr_ = oo_.dr;
|
||||
fname_ = M_.fname;
|
||||
nshock = estim_params_.nvx;
|
||||
nshock = nshock + estim_params_.nvn;
|
||||
nshock = nshock + estim_params_.ncx;
|
||||
nshock = nshock + estim_params_.ncn;
|
||||
|
||||
c0=corrcoef(x);
|
||||
c00=tril(c0,-1);
|
||||
if istab,
|
||||
fig_nam_=[fname_,'_stab_corr_'];
|
||||
else
|
||||
fig_nam_=[fname_,'_unstab_corr_'];
|
||||
end
|
||||
fig_nam_=[fname_,'_',fnam,'_corr_'];
|
||||
|
||||
ifig=0;
|
||||
j2=0;
|
||||
npar=estim_params_.np;
|
||||
if ishock==0
|
||||
npar=estim_params_.np;
|
||||
else
|
||||
npar=estim_params_.np+nshock;
|
||||
end
|
||||
for j=1:npar,
|
||||
i2=find(abs(c00(:,j))>alpha2);
|
||||
if length(i2)>0,
|
||||
|
|
@ -32,11 +45,7 @@ for j=1:npar,
|
|||
j2=j2+1;
|
||||
if mod(j2,12)==1,
|
||||
ifig=ifig+1;
|
||||
if istab
|
||||
figure('name',['Correlations in the stable sample ', num2str(ifig)]),
|
||||
else
|
||||
figure('name',['Correlations in the unstable sample ', num2str(ifig)]),
|
||||
end
|
||||
figure('name',['Correlations in the ',fnam,' sample ', num2str(ifig)]),
|
||||
end
|
||||
subplot(3,4,j2-(ifig-1)*12)
|
||||
% bar(c0(i2,j)),
|
||||
|
|
@ -45,8 +54,15 @@ for j=1:npar,
|
|||
%plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k')
|
||||
%hold on,
|
||||
plot(x(:,j),x(:,i2(jx)),'.')
|
||||
xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'),
|
||||
ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'),
|
||||
% xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'),
|
||||
% ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'),
|
||||
if ishock,
|
||||
xlabel(bayestopt_.name{j},'interpreter','none'),
|
||||
ylabel(bayestopt_.name{i2(jx)},'interpreter','none'),
|
||||
else
|
||||
xlabel(bayestopt_.name{j+nshock},'interpreter','none'),
|
||||
ylabel(bayestopt_.name{i2(jx)+nshock},'interpreter','none'),
|
||||
end
|
||||
title(['cc = ',num2str(c0(i2(jx),j))])
|
||||
if (mod(j2,12)==0) & j2>0,
|
||||
saveas(gcf,[fig_nam_,int2str(ifig)])
|
||||
|
|
@ -58,4 +74,7 @@ for j=1:npar,
|
|||
end
|
||||
|
||||
end
|
||||
if ifig==0,
|
||||
disp(['No correlation term >', num2str(alpha2),' found for ',fnam])
|
||||
end
|
||||
%close all
|
||||
|
|
|
|||
Loading…
Reference in New Issue