diff --git a/matlab/EvalWelfare.m b/matlab/EvalWelfare.m
index dc57638269ef833abc6dc18c22b0c0baf6c182d9..58ff19fe2c9f9634e96aa52cf9cb3e6ea796a5fd 100644
--- a/matlab/EvalWelfare.m
+++ b/matlab/EvalWelfare.m
@@ -46,6 +46,7 @@ windx = strmatch(vName,lgy_(kk,:),'exact');
deep = get_posterior_parameters('posterior_mean');
% [2.3] Compute the posterior distribution of Welfare:
hfid = waitbar(0,'Posterior welfare distribution...');
+compt = 0;
for i=1:B
linea = 1+floor(rand*TotalNumberOfDraws);
tmp = find(METRO(:,3)<linea);
@@ -60,34 +61,36 @@ for i=1:B
eval(['load ' file.name ';'])
DEEP = x2(linee,:);
deep(subindx) = DEEP(subindx);
- deep
set_parameters(deep);
[dr,info] = resol(ys_,0);
if ~info(1)
- WelfDistribution = dr.ys(kk(windx))+.5*dr.ghs2(windx);
+ WelfDistribution(i) = dr.ys(kk(windx))+.5*dr.ghs2(windx);
else
- WelDistribution = Inf;
+ WelfDistribution(i) = Inf;
+ compt = compt+1;
end
waitbar(i/B,hfid);
end
close(hfid);
indx = find(~isinf(WelfDistribution));
-WelfDistribution = WelfDistribution(indx);
-oo_.Welfare.PosteriorDraws = WelfDistribution;
-oo_.Welfare.PosteriorMean = mean(WelfDistribution);
-oo_.Welfare.posteriorStd = std(WelfDistribution);
+WelfareDistribution = WelfDistribution(indx);
+oo_.Welfare.PosteriorDraws = WelfareDistribution;
+oo_.Welfare.PosteriorMean = mean(WelfareDistribution);
+oo_.Welfare.posteriorStd = std(WelfareDistribution);
% [3] Non parametric estimation of the posterior density.
% [3.1] Set some parameters:
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.
% [3.2] Compute the optimal bandwidth parameter:
-optimal_bandwidth = mh_optimal_bandwidth(WelfDistribution,length(indx),bandwidth,kernel_function);
+optimal_bandwidth = mh_optimal_bandwidth(WelfareDistribution,length(indx),bandwidth,kernel_function);
% [3.3] Estimation of the posterior density:
-[x1,f1] = kernel_density_estimate(WelfDistribution,number_of_grid_points,...
+[x1,f1] = kernel_density_estimate(WelfareDistribution,number_of_grid_points,...
optimal_bandwidth,kernel_function);
oo_.Welfare.PosteriorDensity.abscissa = x1;
oo_.Welfare.PosteriorDensity.density = f1;
+disp(['Percentage of mh-draws violating the B&K conditions: ' num2str(compt/B*100) '%.'])
+disp(' ')
% [4] Plot.
if pltOpt
figure('Name','Posterior distribution of the Welfare.')
diff --git a/matlab/metropolis.m b/matlab/metropolis.m
index 94db74a1e1d4a899167608bee68ce096181aa061..9280822ffc61d17d1122f3fe58c4ca1d61a5e807 100644
--- a/matlab/metropolis.m
+++ b/matlab/metropolis.m
@@ -50,9 +50,9 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
if options_.load_mh_file == 0
% Delete old mh files...
if nblck > 1
- disp('MH: Multiple chains mode.')
+ disp('MH: Multiple chains mode.')
else
- disp('MH: One Chain mode.')
+ disp('MH: One Chain mode.')
end
files = eval(['dir(''' fname_ '_mh*.mat'');']);
if length(files)
@@ -115,11 +115,16 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
save([fname_ '_MhInitialization'],'ix2','ilogpo2');
elseif options_.load_mh_file == 1
disp('MH: I''m loading past metropolis-hastings simulations...')
- files = eval(['dir(''' fname_ '_mh*_blck*.mat'');']);
+ if nblck>1
+ files = eval(['dir(''' fname_ '_mh*_blck*.mat'');']);
+ bfiles = eval(['dir(''' fname_ '_mh0_blck*.mat'');']);
+ else
+ files = eval(['dir(''' fname_ '_mh*.mat'');']);
+ bfiles = 1;
+ end
if ~length(files)
error('MH: FAILURE :: there is no MH file to load here!')
end
- bfiles = eval(['dir(''' fname_ '_mh0_blck*.mat'');']);
past_number_of_blocks = length(bfiles);
if length(bfiles)>0 & past_number_of_blocks ~= nblck
disp('MH: The specified number of blocks doesn''t match with the previous number of blocks!')
diff --git a/matlab/subset.m b/matlab/subset.m
index c683092a81935a738e174efdf64cf1f264dbb398..76390fce219f15de182fd152e916d90236873f85 100644
--- a/matlab/subset.m
+++ b/matlab/subset.m
@@ -1,4 +1,4 @@
-function indx = subset();
+function jndx = subset();
% stephane.adjemian@ens.fr [11-30-2005]
global options_ estim_params_ lgx_
@@ -31,22 +31,25 @@ elseif strcmpi(info,'AllWithoutMeasurementErros')
indx = [(1:nvx),nvx+nvn+1:nvx+nvn+ncx,nvx+nvn+ncx+ncn+1:nx]';
end
-
if ~isempty(ExcludedParamNames)
tt = [];
for i = 1:length(ExcludedParamNames)
tmp = strmatch(ExcludedParamNames{i},lgx_);
- if ~isempty(tmp)% The parameter the user wants to exclude is related to the size of the structural innovations.
+ if ~isempty(tmp) & ( strcmpi(info,'All') | strcmpi(info,'StructuralShocks') | ...
+ strcmpi(info,'StructuralShocksWithoutCorrelations') | ...
+ strcmpi(info,'AllWithoutMeasurementErros') )
+ % The parameter the user wants to exclude is related to the size of the structural innovations.
if ncx
disp(['I cannot exclude some the structural variances if the'])
disp(['structural innovations are correlated...'])
error
end
tt = [tt;tmp];
- elseif isempty(tmp) & nvn
+ elseif isempty(tmp) & nvn
tmp = strmatch(ExcludedParamNames{i},options_.varobs);
- if ~isempty(tmp)% The parameter the user wants to exclude is related to the size of the measurement errors.
- % variances:
+ if ~isempty(tmp) & ( strcmpi(info,'All') | strcmpi(info,'MeasurementErrors') | ...
+ strcmpi(info,'MeasurementErrorsWithoutCorrelations') )
+ % The parameter the user wants to exclude is related to the size of the measurement errors variances.
tmp = nvx+tmp;
if ncn
disp(['I cannot exclude some the measurement error variances if the'])
@@ -55,17 +58,20 @@ if ~isempty(ExcludedParamNames)
end
tt = [tt;tmp];
end
- else% Excluded parameters are deep parameters...
- tmp = strmatch(ExcludedParamNames{i},estim_params_.param_names);
+ else% Excluded parameters are deep parameters...
+ tmp = strmatch(ExcludedParamNames{i},estim_params_.param_names,'exact');
if ~isempty(tmp)
- tt = [tt,nvx+nvn+ncx+ncn+tmp];
+ tt = [tt;nvx+nvn+ncx+ncn+tmp];
else
disp('The parameter you want to exclude is unknown!')
error
end
end
end
- for i=1:length(ExcludedParamNames)
- indx = indx(indx~=tt)
+ jndx = [];
+ for i=1:length(indx)
+ if ~any(indx(i)==tt)
+ jndx = [ jndx ; indx(i) ];
+ end
end
end
\ No newline at end of file