diff --git a/matlab/EvalWelfare.m b/matlab/EvalWelfare.m
new file mode 100644
index 0000000000000000000000000000000000000000..dc57638269ef833abc6dc18c22b0c0baf6c182d9
--- /dev/null
+++ b/matlab/EvalWelfare.m
@@ -0,0 +1,99 @@
+function EvalWelfare(vName,pltOpt,NumberOfDraws)
+% stephane.adjemian@ens.fr
+%
+% vName is a string designing the name of the welfare variable in the mod file.
+%
+global options_ oo_ lgy_ ys_ fname_ dr_ info
+
+options_.order = 2;
+options_ = set_default_option(options_,'dr_algo',0);
+options_ = set_default_option(options_,'simul_algo',0);
+subindx = subset();
+
+if isempty(NumberOfDraws)
+ B = 500;
+else
+ B = NumberOfDraws;
+end
+
+% [1] How many mh files and simulations ?
+if options_.mh_nblck == 1
+ files = eval(['dir(''' fname_ '_mh*.mat'');']);
+else% More than one block
+ files = eval(['dir(''' fname_ '_mh*_blck1.mat'');']);
+end
+nn = size(files,1);
+METRO = zeros(nn,3);
+for i = 1:nn
+ file = files(i);
+ eval(['load ' file.name ';'])
+ METRO(i,1) = i-1;
+ METRO(i,2) = size(x2,1);
+ if i>1
+ METRO(i,3) = METRO(i,2)+METRO(i-1,3);
+ else
+ METRO(i,3) = METRO(i,2);
+ end
+end
+TotalNumberOfDraws = METRO(end,3);
+% [2] Take some draws from the metropolis and compute the welfare.
+% [2.1] Set some parameters:
+B = min(B,TotalNumberOfDraws);
+WelfDistribution = zeros(B,1);
+kk = dr_.order_var;
+windx = strmatch(vName,lgy_(kk,:),'exact');
+% [2.2] Get the posterior mean:
+deep = get_posterior_parameters('posterior_mean');
+% [2.3] Compute the posterior distribution of Welfare:
+hfid = waitbar(0,'Posterior welfare distribution...');
+for i=1:B
+ linea = 1+floor(rand*TotalNumberOfDraws);
+ tmp = find(METRO(:,3)<linea);
+ if isempty(tmp)
+ FileNumber = 1;
+ linee = linea;
+ else
+ FileNumber = tmp(end)+1;
+ linee = linea-METRO(tmp(end),3);
+ end
+ file = files(FileNumber);
+ 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);
+ else
+ WelDistribution = Inf;
+ 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);
+% [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);
+% [3.3] Estimation of the posterior density:
+[x1,f1] = kernel_density_estimate(WelfDistribution,number_of_grid_points,...
+ optimal_bandwidth,kernel_function);
+oo_.Welfare.PosteriorDensity.abscissa = x1;
+oo_.Welfare.PosteriorDensity.density = f1;
+% [4] Plot.
+if pltOpt
+ figure('Name','Posterior distribution of the Welfare.')
+ plot(x1,f1,'-k','linewidth',2)
+ axis tight
+ box on
+ xlabel('Welfare')
+ ylabel('Posterior density')
+end
\ No newline at end of file
diff --git a/matlab/metropolis.m b/matlab/metropolis.m
index 1c412edc72f982811892207754b2e129387d81ef..94db74a1e1d4a899167608bee68ce096181aa061 100644
--- a/matlab/metropolis.m
+++ b/matlab/metropolis.m
@@ -1624,11 +1624,10 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
end
if options_.forecast
% FIFTH, I update variable dr_
- resol(ys_,0);
+ dr_ = resol(ys_,0);
% SIXTH, I do and save the forecasts (for all the endogenous variables)
% The state of the economy at the end of the sample
- % depends on the structural parameters.
-
+ % depends on the structural parameters.
yyyy(:,1:ykmin_) = atT(1:endo_nbr,size(atT,2)-ykmin_+1:size(atT,2));
yf = forcst2a(yyyy,dr_,ex_);
if options_.prefilter == 1
@@ -1790,9 +1789,9 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
end
end
if options_.forecast
- dr = resol(ys_,0);
+ dr_ = resol(ys_,0);
yyyy(:,1:ykmin_) = atT(1:endo_nbr,size(atT,2)-ykmin_+1:size(atT,2));
- yf = forcst2a(yyyy,dr,ex_);
+ yf = forcst2a(yyyy,dr_,ex_);
if options_.prefilter == 1
yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
horizon+ykmin_,1);
@@ -1805,7 +1804,7 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
yf = yf+repmat(ys',horizon+ykmin_,1);
end
stock_forcst(:,:,irun_forc) = yf;
- yf1 = forcst2(yyyy,horizon,dr,1);
+ yf1 = forcst2(yyyy,horizon,dr_,1);
if options_.prefilter == 1
yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
repmat(bayestopt_.mean_varobs',[horizon+ykmin_,1,1]);
@@ -1982,7 +1981,7 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
end
end
if options_.forecast
- resol(ys_,0);
+ dr_ = resol(ys_,0);
for j = 1:nvar
yyyy(j,1:ykmin_) = atT(j,size(atT,2)-ykmin_+1:size(atT,2));
end
@@ -2153,79 +2152,79 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
end
end
if options_.forecast
- dr = resol(ys_,0);
- yyyy(:,1:ykmin_) = atT(1:endo_nbr,size(atT,2)-ykmin_+1:size(atT,2));
- yf = forcst2a(yyyy,dr,ex_);
- if options_.prefilter == 1
- yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
- horizon+ykmin_,1);
- end
- yf(:,IdObs) = yf(:,IdObs)+(gend+[1-ykmin_:horizon]')*trend_coeff';
- if options_.loglinear == 1
- yf = yf+repmat(log(ys'),horizon+ykmin_,1);
- yf = exp(yf);
- else
- yf = yf+repmat(ys',horizon+ykmin_,1);
- end
- stock_forcst(:,:,irun_forc) = yf;
- yf1 = forcst2(yyyy,horizon,dr,1);
- if options_.prefilter == 1
- yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
- repmat(bayestopt_.mean_varobs',[horizon+ykmin_,1,1]);
- end
- yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-ykmin_:horizon]')* ...
- trend_coeff',[1,1,1]);
- if options_.loglinear == 1
- yf1 = yf1 + repmat(log(ys'),[horizon+ykmin_,1,1]);
- yf1 = exp(yf1);
- else
- yf1 = yf1 + repmat(ys',[horizon+ykmin_,1,1]);
- end
- stock_forcst1(:,:,irun_forc) = yf1;
- if irun_forc == MAX_nforc
- sfil_forc = sfil_forc + 1;
- save([fname_ '_forecast' int2str(sfil_forc)],'stock_forcst','stock_forcst1');
- irun_forc = 0;
- stock_forcst = zeros(horizon+ykmin_,nvar,MAX_nforc);
- stock_forcst1 = zeros(horizon+ykmin_,nvar,MAX_nforc);
- end
- end
- waitbar(b/B,h);
- end
- if options_.smoother
- if irun_smoo
- stock_smooth = stock_smooth(:,:,1:irun_smoo);
- sfil_smoo = sfil_smoo + 1;
- instr = [fname_ '_smooth' int2str(sfil_smoo) ' stock_smooth;'];
- eval(['save ' instr]);
- end
- clear stock_smooth;
- if irun_inno
- stock_innov = stock_innov(:,:,1:irun_inno);
- sfil_inno = sfil_inno + 1;
- instr = [fname_ '_innovation' int2str(sfil_inno) ' stock_innov;'];
- eval(['save ' instr]);
- end
- clear stock_innov;
- end
- if options_.forecast
- if irun_forc
- stock_forcst = stock_forcst(:,:,1:irun_forc);
- stock_forcst1 = stock_forcst1(:,:,1:irun_forc);
- sfil_forc = sfil_forc + 1;
- save([fname_ '_forecast' int2str(sfil_forc)],'stock_forcst','stock_forcst1');
- end
- clear stock_forcst stock_forcst1
- end
- if options_.filtered_vars
- if irun_filt
- stock_filter = stock_filter(:,:,1:irun_filt);
- sfil_filt = sfil_filt + 1;
- instr = [fname_ '_filter' int2str(sfil_filt) ' stock_filter;'];
- eval(['save ' instr]);
- end
- clear stock_filter;
- end
+ dr_ = resol(ys_,0);
+ yyyy(:,1:ykmin_) = atT(1:endo_nbr,size(atT,2)-ykmin_+1:size(atT,2));
+ yf = forcst2a(yyyy,dr_,ex_);
+ if options_.prefilter == 1
+ yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
+ horizon+ykmin_,1);
+ end
+ yf(:,IdObs) = yf(:,IdObs)+(gend+[1-ykmin_:horizon]')*trend_coeff';
+ if options_.loglinear == 1
+ yf = yf+repmat(log(ys'),horizon+ykmin_,1);
+ yf = exp(yf);
+ else
+ yf = yf+repmat(ys',horizon+ykmin_,1);
+ end
+ stock_forcst(:,:,irun_forc) = yf;
+ yf1 = forcst2(yyyy,horizon,dr_,1);
+ if options_.prefilter == 1
+ yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
+ repmat(bayestopt_.mean_varobs',[horizon+ykmin_,1,1]);
+ end
+ yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-ykmin_:horizon]')* ...
+ trend_coeff',[1,1,1]);
+ if options_.loglinear == 1
+ yf1 = yf1 + repmat(log(ys'),[horizon+ykmin_,1,1]);
+ yf1 = exp(yf1);
+ else
+ yf1 = yf1 + repmat(ys',[horizon+ykmin_,1,1]);
+ end
+ stock_forcst1(:,:,irun_forc) = yf1;
+ if irun_forc == MAX_nforc
+ sfil_forc = sfil_forc + 1;
+ save([fname_ '_forecast' int2str(sfil_forc)],'stock_forcst','stock_forcst1');
+ irun_forc = 0;
+ stock_forcst = zeros(horizon+ykmin_,nvar,MAX_nforc);
+ stock_forcst1 = zeros(horizon+ykmin_,nvar,MAX_nforc);
+ end
+ end
+ waitbar(b/B,h);
+ end
+ if options_.smoother
+ if irun_smoo
+ stock_smooth = stock_smooth(:,:,1:irun_smoo);
+ sfil_smoo = sfil_smoo + 1;
+ instr = [fname_ '_smooth' int2str(sfil_smoo) ' stock_smooth;'];
+ eval(['save ' instr]);
+ end
+ clear stock_smooth;
+ if irun_inno
+ stock_innov = stock_innov(:,:,1:irun_inno);
+ sfil_inno = sfil_inno + 1;
+ instr = [fname_ '_innovation' int2str(sfil_inno) ' stock_innov;'];
+ eval(['save ' instr]);
+ end
+ clear stock_innov;
+ end
+ if options_.forecast
+ if irun_forc
+ stock_forcst = stock_forcst(:,:,1:irun_forc);
+ stock_forcst1 = stock_forcst1(:,:,1:irun_forc);
+ sfil_forc = sfil_forc + 1;
+ save([fname_ '_forecast' int2str(sfil_forc)],'stock_forcst','stock_forcst1');
+ end
+ clear stock_forcst stock_forcst1
+ end
+ if options_.filtered_vars
+ if irun_filt
+ stock_filter = stock_filter(:,:,1:irun_filt);
+ sfil_filt = sfil_filt + 1;
+ instr = [fname_ '_filter' int2str(sfil_filt) ' stock_filter;'];
+ eval(['save ' instr]);
+ end
+ clear stock_filter;
+ end
end
end
close(h);
@@ -2524,7 +2523,7 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
deep(subindx) = DEEP(subindx);
% dsge
set_parameters(deep);
- resol(ys_,0);
+ dr_ = resol(ys_,0);
SS(lgx_orig_ord_,lgx_orig_ord_)=Sigma_e_+1e-14* ...
eye(exo_nbr);
cs = transpose(chol(SS));
@@ -2654,7 +2653,7 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
deep(subindx) = DEEP(subindx);
% dsge
set_parameters(deep);
- resol(ys_,0);
+ dr_ = resol(ys_,0);
SS(lgx_orig_ord_,lgx_orig_ord_) = Sigma_e_+1e-14*eye(exo_nbr);
SS = transpose(chol(SS));
tit(lgx_orig_ord_,:) = lgx_;
@@ -3208,7 +3207,7 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
DEEP = x2(floor(rand*FLN(find(mh_file_number == FLN(:,1)),2))+1,:);
deep(subindx) = DEEP(subindx);
set_parameters(deep);
- resol(ys_,0);
+ dr_ = resol(ys_,0);
Gamma_y = th_autocovariances(dr_,ivar);
if options_.order == 2
m_mean = dr_.ys(ivar) + Gamma_y{options_.ar+3};
@@ -3313,7 +3312,7 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
DEEP = x2(floor(rand*NumberOfSimulations)+1,:);
deep(subindx) = DEEP(subindx);
set_parameters(deep);
- resol(ys_,0);
+ dr_ = resol(ys_,0);
Gamma_y = th_autocovariances(dr_,ivar1);
if options_.order == 2
m_mean = dr_.ys(ivar) + Gamma_y{options_.ar+3};
diff --git a/matlab/resol.m b/matlab/resol.m
index c6e9afd99f6a2e92c583df460f3b2f8150a478be..b53fa179877725546812011c57f65c3d75611213 100644
--- a/matlab/resol.m
+++ b/matlab/resol.m
@@ -5,9 +5,8 @@ function [dr,info]=resol(ys,check_flag)
% plus:
% 11 .... same as dr1 for dr_algo = 2
% 20: can't find steady state info(2) contains sum of sqare residuals
-
- global jacobia_ iy_ ykmin_ ykmax_ gstep_ exo_nbr exo_det_nbr endo_nbr
+global jacobia_ iy_ ykmin_ ykmax_ gstep_ exo_nbr exo_det_nbr endo_nbr
global ex_ ex_det_ valf_ it_ exe_ exe_det_ xkmin_ xkmax_
global fname_ means_ stderrs_ lgy_ maxit_
global dynatol_ options_
diff --git a/matlab/subset.m b/matlab/subset.m
index 3649f369f278e2eca28afcc3e8951fc3282467d5..c683092a81935a738e174efdf64cf1f264dbb398 100644
--- a/matlab/subset.m
+++ b/matlab/subset.m
@@ -17,15 +17,15 @@ nx = nvx+nvn+ncx+ncn+np;
if strcmpi(info,'All')
indx = (1:nx)';
-elseif stcmpi('DeepParameters')
+elseif strcmpi(info,'DeepParameters')
indx = (nvx+nvn+ncx+ncn+1:nx)';
-elseif strcmpi('StructuralShocks')
+elseif strcmpi(info,'StructuralShocks')
indx = [(1:nvx),nvx+nvn+1:nvx+nvn+ncx]';
-elseif stcmpi('StructuralShocksWithoutCorrelations')
+elseif strcmpi(info,'StructuralShocksWithoutCorrelations')
indx = (1:nvx)';
-elseif strcmpi('MeasurementErrors')
+elseif strcmpi(info,'MeasurementErrors')
indx = [(nvx+1:nvx+nvn),(nvx+nvn+ncx+1:nvx+nvn+ncx+ncn)]';
-elseif strcmpi('MeasurementErrorsWithoutCorrelations')
+elseif strcmpi(info,'MeasurementErrorsWithoutCorrelations')
indx = (nvx+1:nvx+nvn)';
elseif strcmpi(info,'AllWithoutMeasurementErros')
indx = [(1:nvx),nvx+nvn+1:nvx+nvn+ncx,nvx+nvn+ncx+ncn+1:nx]';
@@ -34,8 +34,8 @@ end
if ~isempty(ExcludedParamNames)
tt = [];
- for i = 1:length(ParamNames)
- tmp = strmatch(deblank(ExcludedParamNames(i,:)),lgx_);
+ 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 ncx
disp(['I cannot exclude some the structural variances if the'])
@@ -44,7 +44,7 @@ if ~isempty(ExcludedParamNames)
end
tt = [tt;tmp];
elseif isempty(tmp) & nvn
- tmp = strmatch(deblank(ExcludedParamNames(i,:)),options_.varobs);
+ 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:
tmp = nvx+tmp;
@@ -56,7 +56,7 @@ if ~isempty(ExcludedParamNames)
tt = [tt;tmp];
end
else% Excluded parameters are deep parameters...
- tmp = strmatch(deblank(ExcludedParamNames(i,:)),estim_params_.param_names);
+ tmp = strmatch(ExcludedParamNames{i},estim_params_.param_names);
if ~isempty(tmp)
tt = [tt,nvx+nvn+ncx+ncn+tmp];
else
@@ -65,7 +65,7 @@ if ~isempty(ExcludedParamNames)
end
end
end
- for i=1:size(ParamNames,1)
+ for i=1:length(ExcludedParamNames)
indx = indx(indx~=tt)
end
end
\ No newline at end of file