Removed three unused routines.
Stéphane Adjemian (Charybdis)
parent
73cdaa3091
commit
198a6ed7cf
|
|
@ -1,275 +0,0 @@
|
|||
function PosteriorFilterSmootherAndForecast(Y,gend, type,data_index)
|
||||
|
||||
% function PosteriorFilterSmootherAndForecast(Y,gend, type)
|
||||
% Computes posterior filter smoother and forecasts
|
||||
%
|
||||
% INPUTS
|
||||
% Y: data
|
||||
% gend: number of observations
|
||||
% type: posterior
|
||||
% prior
|
||||
% gsa
|
||||
%
|
||||
% OUTPUTS
|
||||
% none
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
|
||||
% Copyright (C) 2005-2012 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
% Dynare is free software: you can redistribute it and/or modify
|
||||
% it under the terms of the GNU General Public License as published by
|
||||
% the Free Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% Dynare is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
global options_ estim_params_ oo_ M_ bayestopt_
|
||||
|
||||
nvx = estim_params_.nvx;
|
||||
nvn = estim_params_.nvn;
|
||||
ncx = estim_params_.ncx;
|
||||
ncn = estim_params_.ncn;
|
||||
np = estim_params_.np ;
|
||||
npar = nvx+nvn+ncx+ncn+np;
|
||||
offset = npar-np;
|
||||
naK = length(options_.filter_step_ahead);
|
||||
|
||||
MaxNumberOfPlotPerFigure = 4;% The square root must be an integer!
|
||||
MaxNumberOfBytes=options_.MaxNumberOfBytes;
|
||||
endo_nbr=M_.endo_nbr;
|
||||
exo_nbr=M_.exo_nbr;
|
||||
nvobs = length(options_.varobs);
|
||||
nn = sqrt(MaxNumberOfPlotPerFigure);
|
||||
iendo = 1:endo_nbr;
|
||||
i_last_obs = gend+(1-M_.maximum_endo_lag:0);
|
||||
horizon = options_.forecast;
|
||||
maxlag = M_.maximum_endo_lag;
|
||||
|
||||
CheckPath('Plots/',M_.dname);
|
||||
MetropolisFolder = CheckPath('metropolis',M_.dname);
|
||||
ModelName = M_.fname;
|
||||
BaseName = [MetropolisFolder filesep Model_name];
|
||||
|
||||
load_last_mh_history_file(MetropolisFolder,ModelName);
|
||||
|
||||
FirstMhFile = record.KeepedDraws.FirstMhFile;
|
||||
FirstLine = record.KeepedDraws.FirstLine;
|
||||
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles;
|
||||
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
|
||||
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
|
||||
|
||||
B = min(1200, round(0.25*NumberOfDraws));
|
||||
B = 200;
|
||||
|
||||
MAX_nruns = min(B,ceil(options_.MaxNumberOfBytes/(npar+2)/8));
|
||||
MAX_nsmoo = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8));
|
||||
MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8));
|
||||
MAX_nerro = min(B,ceil(MaxNumberOfBytes/(length(options_.varobs)*gend)/8));
|
||||
if naK
|
||||
MAX_naK = min(B,ceil(MaxNumberOfBytes/(length(options_.varobs)* ...
|
||||
length(options_.filter_step_ahead)*gend)/8));
|
||||
end
|
||||
if horizon
|
||||
MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8));
|
||||
MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ...
|
||||
8));
|
||||
IdObs = bayestopt_.mfys;
|
||||
|
||||
end
|
||||
|
||||
varlist = options_.varlist;
|
||||
if isempty(varlist)
|
||||
varlist = M_.endo_names;
|
||||
SelecVariables = transpose(1:M_.endo_nbr);
|
||||
nvar = M_.endo_nbr;
|
||||
else
|
||||
nvar = size(varlist,1);
|
||||
SelecVariables = [];
|
||||
for i=1:nvar
|
||||
if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
|
||||
SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
irun1 = 1;
|
||||
irun2 = 1;
|
||||
irun3 = 1;
|
||||
irun4 = 1;
|
||||
irun5 = 1;
|
||||
irun6 = 1;
|
||||
irun7 = 1;
|
||||
ifil1 = 0;
|
||||
ifil2 = 0;
|
||||
ifil3 = 0;
|
||||
ifil4 = 0;
|
||||
ifil5 = 0;
|
||||
ifil6 = 0;
|
||||
ifil7 = 0;
|
||||
|
||||
h = waitbar(0,'Bayesian smoother...');
|
||||
|
||||
stock_param = zeros(MAX_nruns, npar);
|
||||
stock_logpo = zeros(MAX_nruns,1);
|
||||
stock_ys = zeros(MAX_nruns,endo_nbr);
|
||||
if options_.smoother
|
||||
stock_smooth = zeros(endo_nbr,gend,MAX_nsmoo);
|
||||
stock_innov = zeros(exo_nbr,gend,B);
|
||||
stock_error = zeros(nvobs,gend,MAX_nerro);
|
||||
end
|
||||
if options_.filter_step_ahead
|
||||
stock_filter = zeros(naK,endo_nbr,gend+options_.filter_step_ahead(end),MAX_naK);
|
||||
end
|
||||
if options_.forecast
|
||||
stock_forcst_mean = zeros(endo_nbr,horizon+maxlag,MAX_nforc1);
|
||||
stock_forcst_total = zeros(endo_nbr,horizon+maxlag,MAX_nforc2);
|
||||
end
|
||||
|
||||
for b=1:B
|
||||
[deep, logpo] = GetOneDraw(type);
|
||||
M_ = set_all_parameters(deep,estim_params_,M_);
|
||||
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
|
||||
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = ...
|
||||
DsgeSmoother(deep,gend,Y,data_index);
|
||||
|
||||
if options_.loglinear
|
||||
stock_smooth(dr.order_var,:,irun1) = alphahat(1:endo_nbr,:)+ ...
|
||||
repmat(log(dr.ys(dr.order_var)),1,gend);
|
||||
else
|
||||
stock_smooth(dr.order_var,:,irun1) = alphahat(1:endo_nbr,:)+ ...
|
||||
repmat(dr.ys(dr.order_var),1,gend);
|
||||
end
|
||||
if nvx
|
||||
stock_innov(:,:,irun2) = etahat;
|
||||
end
|
||||
if nvn
|
||||
stock_error(:,:,irun3) = epsilonhat;
|
||||
end
|
||||
if naK
|
||||
stock_filter(:,dr.order_var,:,irun4) = aK(options_.filter_step_ahead,1:endo_nbr,:);
|
||||
end
|
||||
stock_param(irun5,:) = deep;
|
||||
stock_logpo(irun5,1) = logpo;
|
||||
stock_ys(irun5,:) = SteadyState';
|
||||
|
||||
if horizon
|
||||
yyyy = alphahat(iendo,i_last_obs);
|
||||
yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
|
||||
if options_.prefilter == 1
|
||||
yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
|
||||
horizon+maxlag,1);
|
||||
end
|
||||
yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff';
|
||||
if options_.loglinear
|
||||
yf = yf+repmat(log(SteadyState'),horizon+maxlag,1);
|
||||
else
|
||||
yf = yf+repmat(SteadyState',horizon+maxlag,1);
|
||||
end
|
||||
yf1 = forcst2(yyyy,horizon,dr,1);
|
||||
if options_.prefilter
|
||||
yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
|
||||
repmat(bayestopt_.mean_varobs',[horizon+maxlag,1,1]);
|
||||
end
|
||||
yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-maxlag:horizon]')* ...
|
||||
trend_coeff',[1,1,1]);
|
||||
if options_.loglinear
|
||||
yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
|
||||
else
|
||||
yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]);
|
||||
end
|
||||
|
||||
stock_forcst_mean(:,:,irun6) = yf';
|
||||
stock_forcst_total(:,:,irun7) = yf1';
|
||||
end
|
||||
|
||||
irun1 = irun1 + 1;
|
||||
irun2 = irun2 + 1;
|
||||
irun3 = irun3 + 1;
|
||||
irun4 = irun4 + 1;
|
||||
irun5 = irun5 + 1;
|
||||
irun6 = irun6 + 1;
|
||||
irun7 = irun7 + 1;
|
||||
|
||||
if irun1 > MAX_nsmoo || b == B
|
||||
stock = stock_smooth(:,:,1:irun1-1);
|
||||
ifil1 = ifil1 + 1;
|
||||
save([BaseName '_smooth' int2str(ifil1) '.mat'],'stock');
|
||||
irun1 = 1;
|
||||
end
|
||||
|
||||
if nvx && (irun2 > MAX_ninno || b == B)
|
||||
stock = stock_innov(:,:,1:irun2-1);
|
||||
ifil2 = ifil2 + 1;
|
||||
save([BaseName '_inno' int2str(ifil2) '.mat'],'stock');
|
||||
irun2 = 1;
|
||||
end
|
||||
|
||||
if nvn && (irun3 > MAX_error || b == B)
|
||||
stock = stock_error(:,:,1:irun3-1);
|
||||
ifil3 = ifil3 + 1;
|
||||
save([BaseName '_error' int2str(ifil3) '.mat'],'stock');
|
||||
irun3 = 1;
|
||||
end
|
||||
|
||||
if naK && (irun4 > MAX_naK || b == B)
|
||||
stock = stock_filter(:,:,:,1:irun4-1);
|
||||
ifil4 = ifil4 + 1;
|
||||
save([BaseName '_filter' int2str(ifil4) '.mat'],'stock');
|
||||
irun4 = 1;
|
||||
end
|
||||
|
||||
if irun5 > MAX_nruns || b == B
|
||||
stock = stock_param(1:irun5-1,:);
|
||||
ifil5 = ifil5 + 1;
|
||||
save([BaseName '_param' int2str(ifil5) '.mat'],'stock','stock_logpo','stock_ys');
|
||||
irun5 = 1;
|
||||
end
|
||||
|
||||
if horizon && (irun6 > MAX_nforc1 || b == B)
|
||||
stock = stock_forcst_mean(:,:,1:irun6-1);
|
||||
ifil6 = ifil6 + 1;
|
||||
save([BaseName '_forc_mean' int2str(ifil6) '.mat'],'stock');
|
||||
irun6 = 1;
|
||||
end
|
||||
|
||||
if horizon && (irun7 > MAX_nforc2 || b == B)
|
||||
stock = stock_forcst_total(:,:,1:irun7-1);
|
||||
ifil7 = ifil7 + 1;
|
||||
save([BaseName '_forc_total' int2str(ifil7) '.mat'],'stock');
|
||||
irun7 = 1;
|
||||
end
|
||||
|
||||
waitbar(b/B,h);
|
||||
|
||||
end
|
||||
close(h)
|
||||
|
||||
% ??????
|
||||
stock_gend=gend;
|
||||
stock_data=Y;
|
||||
save([BaseName '_data.mat'],'stock_gend','stock_data');
|
||||
|
||||
if options_.smoother
|
||||
pm3(endo_nbr,gend,ifil1,B,'Smoothed variables',...
|
||||
M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
|
||||
'names2','smooth',MetropolisFolder,'_smooth')
|
||||
end
|
||||
|
||||
if options_.forecast
|
||||
pm3(endo_nbr,horizon+maxlag,ifil6,B,'Forecasted variables (mean)',...
|
||||
M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
|
||||
'names2','smooth',MetropolisFolder,'_forc_mean')
|
||||
pm3(endo_nbr,horizon+maxlag,ifil6,B,'Forecasted variables (total)',...
|
||||
M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
|
||||
'names2','smooth',MetropolisFolder,'_forc_total')
|
||||
end
|
||||
|
|
@ -1,152 +0,0 @@
|
|||
function forcst_unc(y0,var_list)
|
||||
% function [mean,intval1,intval2]=forcst_unc(y0,var_list)
|
||||
% computes forecasts with parameter uncertainty
|
||||
%
|
||||
% INPUTS
|
||||
% y0: matrix of initial values
|
||||
% var_list: list of variables to be forecasted
|
||||
%
|
||||
% OUTPUTS
|
||||
% none
|
||||
%
|
||||
% ALGORITHM
|
||||
% uses antithetic draws for the shocks
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% None.
|
||||
|
||||
% Copyright (C) 2006-2013 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
% Dynare is free software: you can redistribute it and/or modify
|
||||
% it under the terms of the GNU General Public License as published by
|
||||
% the Free Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% Dynare is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% 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_ options_ oo_ estim_params_ bayestopt_
|
||||
|
||||
% setting up estim_params_
|
||||
[xparam1,estim_params_,bayestopt_,lb,ub] = set_prior(estim_params_,M_);
|
||||
|
||||
options_.TeX = 0;
|
||||
options_.nograph = 0;
|
||||
plot_priors(bayestopt_,M_,options_);
|
||||
|
||||
% workspace initialization
|
||||
if isempty(var_list)
|
||||
var_list = M_.endo_names(1:M_.orig_endo_nbr,:);
|
||||
end
|
||||
n = size(var_list,1);
|
||||
|
||||
periods = options_.forecast;
|
||||
exo_nbr = M_.exo_nbr;
|
||||
replic = options_.replic;
|
||||
order = options_.order;
|
||||
maximum_lag = M_.maximum_lag;
|
||||
% params = prior_draw(1);
|
||||
params = rndprior(bayestopt_);
|
||||
set_parameters(params);
|
||||
% eliminate shocks with 0 variance
|
||||
i_exo_var = setdiff([1:exo_nbr],find(diag(M_.Sigma_e) == 0));
|
||||
nx = length(i_exo_var);
|
||||
|
||||
ex0 = zeros(periods,exo_nbr);
|
||||
yf1 = zeros(periods+M_.maximum_lag,n,replic);
|
||||
|
||||
% loops on parameter values
|
||||
m1 = 0;
|
||||
m2 = 0;
|
||||
for i=1:replic
|
||||
% draw parameter values from the prior
|
||||
% params = prior_draw(0);
|
||||
params = rndprior(bayestopt_);
|
||||
set_parameters(params);
|
||||
% solve the model
|
||||
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
|
||||
% discard problematic cases
|
||||
if info
|
||||
continue
|
||||
end
|
||||
% compute forecast with zero shocks
|
||||
m1 = m1+1;
|
||||
yf1(:,:,m1) = simult_(y0,dr,ex0,order)';
|
||||
% compute forecast with antithetic shocks
|
||||
chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var));
|
||||
ex(:,i_exo_var) = randn(periods,nx)*chol_S;
|
||||
m2 = m2+1;
|
||||
yf2(:,:,m2) = simult_(y0,dr,ex,order)';
|
||||
m2 = m2+1;
|
||||
yf2(:,:,m2) = simult_(y0,dr,-ex,order)';
|
||||
end
|
||||
|
||||
oo_.forecast.accept_rate = (replic-m1)/replic;
|
||||
|
||||
if options_.noprint == 0 && m1 < replic
|
||||
skipline(2)
|
||||
disp('FORECASTING WITH PARAMETER UNCERTAINTY:')
|
||||
disp(sprintf(['The model couldn''t be solved for %f%% of the parameter' ...
|
||||
' values'],100*oo_.forecast.accept_rate))
|
||||
skipline(2)
|
||||
end
|
||||
|
||||
% compute moments
|
||||
yf1 = yf1(:,:,1:m1);
|
||||
yf2 = yf2(:,:,1:m2);
|
||||
|
||||
yf_mean = mean(yf1,3);
|
||||
|
||||
yf1 = sort(yf1,3);
|
||||
yf2 = sort(yf2,3);
|
||||
|
||||
sig_lev = options_.conf_sig;
|
||||
k1 = round((0.5+[-sig_lev, sig_lev]/2)*replic);
|
||||
% make sure that lower bound is at least the first element
|
||||
if k1(2) == 0
|
||||
k1(2) = 1;
|
||||
end
|
||||
k2 = round((1+[-sig_lev, sig_lev])*replic);
|
||||
% make sure that lower bound is at least the first element
|
||||
if k2(2) == 0
|
||||
k2(2) = 1;
|
||||
end
|
||||
|
||||
% compute shock uncertainty around forecast with mean prior
|
||||
set_parameters(bayestopt_.p1);
|
||||
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
|
||||
[yf3,yf3_intv] = forcst(dr,y0,periods,var_list);
|
||||
yf3_1 = yf3'-[zeros(maximum_lag,n); yf3_intv];
|
||||
yf3_2 = yf3'+[zeros(maximum_lag,n); yf3_intv];
|
||||
|
||||
% graphs
|
||||
OutputDirectoryName = CheckPath('graphs',M_.fname);
|
||||
|
||||
dyn_graph=dynare_graph_init('Forecasts type I',n,{'b-' 'g-' 'g-' 'r-' 'r-'});
|
||||
for i=1:n
|
||||
dynare_graph(dyn_graph,[yf_mean(:,i) squeeze(yf1(:,i,k1)) squeeze(yf2(:,i,k2))],...
|
||||
var_list(i,:));
|
||||
end
|
||||
dyn_saveas(dyn_graph.fh,[OutputDirectoryName '/' M_.fname '_forecast_param_uncert_',num2str(nlags)],options_)
|
||||
|
||||
dyn_graph=dynare_graph_init('Forecasts type II',n,{'b-' 'k-' 'k-' 'r-' 'r-'});
|
||||
for i=1:n
|
||||
dynare_graph(dyn_graph,[yf_mean(:,i) yf3_1(:,i) yf3_2(:,i) squeeze(yf2(:,i,k2))],...
|
||||
var_list(i,:));
|
||||
end
|
||||
dyn_saveas(dyn_graph.fh,[OutputDirectoryName '/' M_.fname '_forecast_param_shock_uncert_',num2str(nlags)],options_)
|
||||
|
||||
|
||||
% saving results
|
||||
save_results(yf_mean,'oo_.forecast.Mean.',var_list);
|
||||
save_results(yf1(:,:,k1(1)),'oo_.forecast.HPDinf.',var_list);
|
||||
save_results(yf1(:,:,k1(2)),'oo_.forecast.HPDsup.',var_list);
|
||||
save_results(yf2(:,:,k2(1)),'oo_.forecast.HPDTotalinf.',var_list);
|
||||
save_results(yf2(:,:,k2(2)),'oo_.forecast.HPDTotalsup.',var_list);
|
||||
|
|
@ -1,57 +0,0 @@
|
|||
function y = rndprior(bayestopt_)
|
||||
% function y = rndprior(bayestopt_)
|
||||
% Draws random number from the prior density
|
||||
%
|
||||
% INPUTS
|
||||
% bayestopt_: structure characterizing priors
|
||||
%
|
||||
% OUTPUTS
|
||||
% y: drawn numbers vector
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
|
||||
% Copyright (C) 2003-2009 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
% Dynare is free software: you can redistribute it and/or modify
|
||||
% it under the terms of the GNU General Public License as published by
|
||||
% the Free Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% Dynare is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
pshape=bayestopt_.pshape;
|
||||
p3=bayestopt_.p3;
|
||||
p4=bayestopt_.p4;
|
||||
p6=bayestopt_.p6;
|
||||
p7=bayestopt_.p7;
|
||||
|
||||
y = NaN(1,length(pshape));
|
||||
|
||||
for i=1:length(pshape)
|
||||
switch pshape(i)
|
||||
case 1 % Beta
|
||||
y(i) = betarnd(p6(i),p7(i));
|
||||
y(i) = y(1,i) * (p4(i)-p3(i)) + p3(i);
|
||||
case 2 % Generalized gamma
|
||||
y(i) = gamrnd(p6(i),p7(i)) + p3(i);
|
||||
case 3 % Gaussian
|
||||
y(i) = randn*p7(i) + p6(i) ;
|
||||
case 4 % Inverse-gamma type 1
|
||||
y(i) = 1/sqrt(gamrnd(p7(i)/2, 2/p6(i))) + p3(i);
|
||||
case 5 % Uniform
|
||||
y(i) = rand*(p4(i)-p3(i)) + p3(i);
|
||||
case 6 % Inverse-gamma type 2
|
||||
y(i) = 1/gamrnd(p7(i)/2, 2/p6(i)) + p3(i);
|
||||
otherwise
|
||||
error(sprintf('rndprior: unknown distribution shape (index %d, type %d)', i, pshape(i)));
|
||||
end
|
||||
end
|
||||
Loading…
Reference in New Issue