Consider measurement error for uncertainty bands in Bayesian estimation

matlab/forcst2.m

@@ -1,7 +1,8 @@
 function yf=forcst2(y0,horizon,dr,n)
 % function yf=forcst2(y0,horizon,dr,n)
 % 
-% computes forecasts based on first order model solution, given shocks drawn from the shock distribution
+% computes forecasts based on first order model solution, given shocks
+% drawn from the shock distribution, but not including measurement error
 % Inputs:
 %   - y0        [endo_nbr by maximum_endo_lag]          matrix of starting values
 %   - dr        [structure]                         structure with Dynare decision rules
@@ -11,7 +12,7 @@ function yf=forcst2(y0,horizon,dr,n)
 % Outputs:
 %   - yf        [horizon+ykmin_ by endo_nbr by n]   array of forecasts
 % 
-% Copyright (C) 2008-2015 Dynare Team
+% Copyright (C) 2008-2016 Dynare Team
 %
 % This file is part of Dynare.
 %

matlab/prior_posterior_statistics.m

@@ -107,6 +107,9 @@ 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));
+    if ~isequal(M_.H,0)
+        MAX_nforc_ME = min(B,ceil(MaxNumberOfBytes/((size(options_.varobs,1))*(horizon+maxlag))/8));
+    end
 end
 MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8)));
 
@@ -126,7 +129,7 @@ for i=1:nvar
     end
 end
 
-n_variables_to_fill=11;
+n_variables_to_fill=12;
 
 irun = ones(n_variables_to_fill,1);
 ifil = zeros(n_variables_to_fill,1);
@@ -163,6 +166,9 @@ if horizon
     localVars.i_last_obs=i_last_obs;
     localVars.MAX_nforc1=MAX_nforc1;
     localVars.MAX_nforc2=MAX_nforc2;
+    if ~isequal(M_.H,0)
+        localVars.MAX_nforc_ME = MAX_nforc_ME;
+    end
 end
 localVars.exo_nbr=exo_nbr;
 localVars.maxlag=maxlag;
@@ -224,6 +230,10 @@ else
             ifil(6,j+1) =ifil(6,j)+nfiles;
             nfiles = ceil(nBlockPerCPU(j)/MAX_nforc2);
             ifil(7,j+1) =ifil(7,j)+nfiles;
+            if ~isequal(M_.H,0)
+                nfiles = ceil(nBlockPerCPU(j)/MAX_nforc_ME);
+                ifil(12,j+1) =ifil(12,j)+nfiles;
+            end
         end
         if options_.filter_covariance
             nfiles = ceil(nBlockPerCPU(j)/MAX_filter_covariance);
@@ -318,6 +328,19 @@ if options_.forecast
     pm3(endo_nbr,horizon,ifil(7),B,'Forecasted variables (point)',...
         '',varlist,M_.endo_names_tex,M_.endo_names,...
         varlist,'PointForecast',DirectoryName,'_forc_point');
+    if ~isequal(M_.H,0)
+        texnames=[];
+        for obs_iter=1:length(options_.varobs)        
+            obs_names{obs_iter,1}=M_.endo_names(strmatch(options_.varobs{obs_iter},M_.endo_names,'exact'),:);
+            texnames{obs_iter,1}=M_.endo_names_tex(strmatch(options_.varobs{obs_iter},M_.endo_names,'exact'),:);
+        end
+        obs_names=char(obs_names);
+        texnames=char(texnames);
+        varlist_forecast_ME=intersect(options_.varobs,varlist);
+        pm3(meas_err_nbr,horizon,ifil(12),B,'Forecasted variables (point) with ME',...
+           '',char(varlist_forecast_ME),texnames,obs_names,...
+           char(varlist_forecast_ME),'PointForecastME',DirectoryName,'_forc_point_ME')    
+    end
 end
 
 if options_.filter_covariance

matlab/prior_posterior_statistics_core.m

@@ -16,6 +16,7 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa
 %                          _param;
 %                          _forc_mean;
 %                          _forc_point;
+%                          _forc_point_ME;
 %                          _filter_covar;
 %                          _trend_coeff;
 %                          _smoothed_trend;
@@ -74,6 +75,9 @@ if horizon
     i_last_obs=myinputs.i_last_obs;
     MAX_nforc1=myinputs.MAX_nforc1;
     MAX_nforc2=myinputs.MAX_nforc2;
+    if ~isequal(M_.H,0)
+        MAX_nforc_ME=myinputs.MAX_nforc_ME;
+    end
 end
 if naK
     MAX_naK=myinputs.MAX_naK;
@@ -137,6 +141,7 @@ if RemoteFlag==1,
     OutputFileName_param = {};
     OutputFileName_forc_mean = {};
     OutputFileName_forc_point = {};
+    OutputFileName_forc_point_ME = {};
     OutputFileName_filter_covar ={};
     OutputFileName_trend_coeff = {};
     OutputFileName_smoothed_trend = {};
@@ -155,6 +160,9 @@ if run_smoother
   if horizon
       stock_forcst_mean= NaN(endo_nbr,horizon,MAX_nforc1);
       stock_forcst_point = NaN(endo_nbr,horizon,MAX_nforc2);
+      if ~isequal(M_.H,0)
+          stock_forcst_point_ME = NaN(length(varobs),horizon,MAX_nforc_ME);
+      end
   end
 end
 if nvn
@@ -300,6 +308,16 @@ for b=fpar:B
 
             stock_forcst_mean(:,:,irun(6)) = yf(maxlag+1:end,:)';
             stock_forcst_point(:,:,irun(7)) = yf1(maxlag+1:end,:)';
+            if ~isequal(M_.H,0)
+                ME_shocks=zeros(length(varobs),horizon);
+                i_exo_var = setdiff([1:length(varobs)],find(diag(M_.H) == 0));
+                nxs = length(i_exo_var);
+                chol_H = chol(M_.H(i_exo_var,i_exo_var));
+                if ~isempty(M_.H)
+                    ME_shocks(i_exo_var,:) = chol_H*randn(nxs,horizon);
+                end
+                stock_forcst_point_ME(:,:,irun(12)) = yf1(maxlag+1:end,IdObs)'+ME_shocks;
+            end
         end
         if filter_covariance
             stock_filter_covariance(dr.order_var,dr.order_var,:,irun(8)) = P;
@@ -312,7 +330,7 @@ for b=fpar:B
     stock_ys(irun(5),:) = SteadyState';
 
 
-    irun = irun +  ones(11,1);
+    irun = irun +  ones(12,1);
 
 
     if run_smoother && (irun(1) > MAX_nsmoo || b == B),
@@ -433,7 +451,18 @@ for b=fpar:B
         end
         irun(irun_index) = 1;
     end
-
+    
+    irun_index=12;
+    if run_smoother && horizon && ~isequal(M_.H,0) && (irun(irun_index) > MAX_nforc_ME ||  b == B)
+        stock = stock_forcst_point_ME(:,:,1:irun(irun_index)-1);
+        ifil(irun_index) = ifil(irun_index) + 1;
+        save([DirectoryName '/' M_.fname '_forc_point_ME' int2str(ifil(irun_index)) '.mat'],'stock');
+        if RemoteFlag==1,
+            OutputFileName_forc_point_ME = [OutputFileName_forc_point_ME; {[DirectoryName filesep], [M_.fname '_forc_point_ME' int2str(ifil(irun_index)) '.mat']}];
+        end
+        irun(irun_index) = 1;
+    end
+    
     dyn_waitbar((b-fpar+1)/(B-fpar+1),h);
 end
 
@@ -447,6 +476,7 @@ if RemoteFlag==1,
                         OutputFileName_param;
                         OutputFileName_forc_mean;
                         OutputFileName_forc_point;
+                        OutputFileName_forc_point_ME;
                         OutputFileName_filter_covar;
                         OutputFileName_trend_coeff;
                         OutputFileName_smoothed_trend;