Description of PosteriorFilterSmootherAndForecast

0001 function PosteriorFilterSmootherAndForecast(Y,gend, type,data_index)
0002 
0003 % function PosteriorFilterSmootherAndForecast(Y,gend, type)
0004 % Computes posterior filter smoother and forecasts
0005 %
0006 % INPUTS
0007 %    Y:      data
0008 %    gend:   number of observations
0009 %    type:   posterior
0010 %            prior
0011 %            gsa
0012 %
0013 % OUTPUTS
0014 %    none
0015 %
0016 % SPECIAL REQUIREMENTS
0017 %    none
0018 
0019 % Copyright (C) 2005-2011 Dynare Team
0020 %
0021 % This file is part of Dynare.
0022 %
0023 % Dynare is free software: you can redistribute it and/or modify
0024 % it under the terms of the GNU General Public License as published by
0025 % the Free Software Foundation, either version 3 of the License, or
0026 % (at your option) any later version.
0027 %
0028 % Dynare is distributed in the hope that it will be useful,
0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0031 % GNU General Public License for more details.
0032 %
0033 % You should have received a copy of the GNU General Public License
0034 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0035 
0036 global options_ estim_params_ oo_ M_ bayestopt_
0037 
0038 nvx  = estim_params_.nvx;
0039 nvn  = estim_params_.nvn;
0040 ncx  = estim_params_.ncx;
0041 ncn  = estim_params_.ncn;
0042 np   = estim_params_.np ;
0043 npar = nvx+nvn+ncx+ncn+np;
0044 offset = npar-np;
0045 naK = length(options_.filter_step_ahead);
0046 %%
0047 MaxNumberOfPlotPerFigure = 4;% The square root must be an integer!
0048 MaxNumberOfBytes=options_.MaxNumberOfBytes;
0049 endo_nbr=M_.endo_nbr;
0050 exo_nbr=M_.exo_nbr;
0051 nvobs     = size(options_.varobs,1);
0052 nn = sqrt(MaxNumberOfPlotPerFigure);
0053 iendo = 1:endo_nbr;
0054 i_last_obs = gend+(1-M_.maximum_endo_lag:0);
0055 horizon = options_.forecast;
0056 maxlag = M_.maximum_endo_lag;
0057 %%
0058 CheckPath('Plots/',M_.dname);
0059 DirectoryName = CheckPath('metropolis',M_.dname);
0060 load([ DirectoryName '/'  M_.fname '_mh_history.mat'])
0061 FirstMhFile = record.KeepedDraws.FirstMhFile;
0062 FirstLine = record.KeepedDraws.FirstLine;
0063 TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles;
0064 TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
0065 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
0066 clear record;
0067 B = min(1200, round(0.25*NumberOfDraws));
0068 B = 200;
0069 %%
0070 MAX_nruns = min(B,ceil(options_.MaxNumberOfBytes/(npar+2)/8));
0071 MAX_nsmoo = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8));
0072 MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8));
0073 MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8));
0074 if naK
0075     MAX_naK   = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)* ...
0076                                              length(options_.filter_step_ahead)*gend)/8));
0077 end
0078 if horizon
0079     MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8));
0080     MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ...
0081                             8));
0082     IdObs    = bayestopt_.mfys;
0083 
0084 end
0085 
0086 %%
0087 varlist = options_.varlist;
0088 if isempty(varlist)
0089     varlist = M_.endo_names;
0090     SelecVariables = transpose(1:M_.endo_nbr);
0091     nvar = M_.endo_nbr;
0092 else
0093     nvar = size(varlist,1);
0094     SelecVariables = [];
0095     for i=1:nvar
0096         if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
0097             SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
0098         end
0099     end
0100 end
0101 
0102 irun1 = 1;
0103 irun2 = 1;
0104 irun3 = 1;
0105 irun4 = 1;
0106 irun5 = 1;
0107 irun6 = 1;
0108 irun7 = 1;
0109 ifil1 = 0;
0110 ifil2 = 0;
0111 ifil3 = 0;
0112 ifil4 = 0;
0113 ifil5 = 0;
0114 ifil6 = 0;
0115 ifil7 = 0;
0116 
0117 h = waitbar(0,'Bayesian smoother...');
0118 
0119 stock_param = zeros(MAX_nruns, npar);
0120 stock_logpo = zeros(MAX_nruns,1);
0121 stock_ys = zeros(MAX_nruns,endo_nbr);
0122 if options_.smoother
0123     stock_smooth = zeros(endo_nbr,gend,MAX_nsmoo);
0124     stock_innov  = zeros(exo_nbr,gend,B);
0125     stock_error = zeros(nvobs,gend,MAX_nerro);
0126 end
0127 if options_.filter_step_ahead
0128     stock_filter = zeros(naK,endo_nbr,gend+options_.filter_step_ahead(end),MAX_naK);
0129 end
0130 if options_.forecast
0131     stock_forcst_mean = zeros(endo_nbr,horizon+maxlag,MAX_nforc1);
0132     stock_forcst_total = zeros(endo_nbr,horizon+maxlag,MAX_nforc2);
0133 end
0134 
0135 for b=1:B
0136     %deep = GetOneDraw(NumberOfDraws,FirstMhFile,LastMhFile,FirstLine,MAX_nruns,DirectoryName);
0137     [deep, logpo] = GetOneDraw(type);
0138     set_all_parameters(deep);
0139     [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
0140     [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = ...
0141         DsgeSmoother(deep,gend,Y,data_index);
0142 
0143     if options_.loglinear
0144         stock_smooth(dr.order_var,:,irun1) = alphahat(1:endo_nbr,:)+ ...
0145             repmat(log(dr.ys(dr.order_var)),1,gend);
0146     else
0147         stock_smooth(dr.order_var,:,irun1) = alphahat(1:endo_nbr,:)+ ...
0148             repmat(dr.ys(dr.order_var),1,gend);
0149     end
0150     if nvx
0151         stock_innov(:,:,irun2)  = etahat;
0152     end
0153     if nvn
0154         stock_error(:,:,irun3)  = epsilonhat;
0155     end
0156     if naK
0157         stock_filter(:,dr.order_var,:,irun4) = aK(options_.filter_step_ahead,1:endo_nbr,:);
0158     end
0159     stock_param(irun5,:) = deep;
0160     stock_logpo(irun5,1) = logpo;
0161     stock_ys(irun5,:) = SteadyState';
0162 
0163     if horizon
0164         yyyy = alphahat(iendo,i_last_obs);
0165         yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
0166         if options_.prefilter == 1
0167             yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
0168                                              horizon+maxlag,1);
0169         end
0170         yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff';
0171         if options_.loglinear == 1
0172             yf = yf+repmat(log(SteadyState'),horizon+maxlag,1);
0173             %      yf = exp(yf);
0174         else
0175             yf = yf+repmat(SteadyState',horizon+maxlag,1);
0176         end
0177         yf1 = forcst2(yyyy,horizon,dr,1);
0178         if options_.prefilter == 1
0179             yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
0180                 repmat(bayestopt_.mean_varobs',[horizon+maxlag,1,1]);
0181         end
0182         yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-maxlag:horizon]')* ...
0183                                                trend_coeff',[1,1,1]);
0184         if options_.loglinear == 1
0185             yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
0186             %      yf1 = exp(yf1);
0187         else
0188             yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]);
0189         end
0190 
0191         stock_forcst_mean(:,:,irun6) = yf';
0192         stock_forcst_total(:,:,irun7) = yf1';
0193     end
0194 
0195     irun1 = irun1 + 1;
0196     irun2 = irun2 + 1;
0197     irun3 = irun3 + 1;
0198     irun4 = irun4 + 1;
0199     irun5 = irun5 + 1;
0200     irun6 = irun6 + 1;
0201     irun7 = irun7 + 1;
0202 
0203     if irun1 > MAX_nsmoo || b == B
0204         stock = stock_smooth(:,:,1:irun1-1);
0205         ifil1 = ifil1 + 1;
0206         save([DirectoryName '/' M_.fname '_smooth' int2str(ifil1) '.mat'],'stock');
0207         irun1 = 1;
0208     end
0209 
0210     if nvx && (irun2 > MAX_ninno || b == B)
0211         stock = stock_innov(:,:,1:irun2-1);
0212         ifil2 = ifil2 + 1;
0213         save([DirectoryName '/' M_.fname '_inno' int2str(ifil2) '.mat'],'stock');
0214         irun2 = 1;
0215     end
0216 
0217     if nvn && (irun3 > MAX_error || b == B)
0218         stock = stock_error(:,:,1:irun3-1);
0219         ifil3 = ifil3 + 1;
0220         save([DirectoryName '/' M_.fname '_error' int2str(ifil3) '.mat'],'stock');
0221         irun3 = 1;
0222     end
0223 
0224     if naK && (irun4 > MAX_naK || b == B)
0225         stock = stock_filter(:,:,:,1:irun4-1);
0226         ifil4 = ifil4 + 1;
0227         save([DirectoryName '/' M_.fname '_filter' int2str(ifil4) '.mat'],'stock');
0228         irun4 = 1;
0229     end
0230 
0231     if irun5 > MAX_nruns || b == B
0232         stock = stock_param(1:irun5-1,:);
0233         ifil5 = ifil5 + 1;
0234         save([DirectoryName '/' M_.fname '_param' int2str(ifil5) '.mat'],'stock','stock_logpo','stock_ys');
0235         irun5 = 1;
0236     end
0237 
0238     if horizon && (irun6 > MAX_nforc1 || b == B)
0239         stock = stock_forcst_mean(:,:,1:irun6-1);
0240         ifil6 = ifil6 + 1;
0241         save([DirectoryName '/' M_.fname '_forc_mean' int2str(ifil6) '.mat'],'stock');
0242         irun6 = 1;
0243     end
0244 
0245     if horizon && (irun7 > MAX_nforc2 ||  b == B)
0246         stock = stock_forcst_total(:,:,1:irun7-1);
0247         ifil7 = ifil7 + 1;
0248         save([DirectoryName '/' M_.fname '_forc_total' int2str(ifil7) '.mat'],'stock');
0249         irun7 = 1;
0250     end
0251 
0252     waitbar(b/B,h);
0253 end
0254 close(h)
0255 
0256 stock_gend=gend;
0257 stock_data=Y;
0258 save([DirectoryName '/' M_.fname '_data.mat'],'stock_gend','stock_data');
0259 
0260 if options_.smoother
0261     pm3(endo_nbr,gend,ifil1,B,'Smoothed variables',...
0262         M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
0263         'names2','smooth',[M_.fname '/metropolis'],'_smooth')
0264 end
0265 
0266 if options_.forecast
0267     pm3(endo_nbr,horizon+maxlag,ifil6,B,'Forecasted variables (mean)',...
0268         M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
0269         'names2','smooth',[M_.fname '/metropolis'],'_forc_mean')
0270     pm3(endo_nbr,horizon+maxlag,ifil6,B,'Forecasted variables (total)',...
0271         M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
0272         'names2','smooth',[M_.fname '/metropolis'],'_forc_total')
0273 end
PosteriorFilterSmootherAndForecast

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
