Fixed bug introduced with the parallelization of PosteriorIRF.m.

Changes related to the new interface for dsge-var models.

matlab/PosteriorIRF.m

@@ -62,48 +62,53 @@ ncx  = estim_params_.ncx;
 ncn  = estim_params_.ncn;
 np   = estim_params_.np ;
 npar = nvx+nvn+ncx+ncn+np;
-offset = npar-np;
-clear('nvx','nvn','ncx','ncn','np');
+offset = npar-np; clear('nvx','nvn','ncx','ncn','np');
+
 nvobs = size(options_.varobs,1);
 gend = options_.nobs;
 MaxNumberOfPlotPerFigure = 9;
 nn = sqrt(MaxNumberOfPlotPerFigure);
 MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8);
 MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
-if ~isempty(strmatch('dsge_prior_weight',M_.param_names))
+if options_.dsge_var
     MAX_nirfs_dsgevar = ceil(options_.MaxNumberOfBytes/(options_.irf*nvobs*M_.exo_nbr)/8);
 else
     MAX_nirfs_dsgevar = 0;
 end
+
 DirectoryName = CheckPath('Output');
 if strcmpi(type,'posterior')
-  MhDirectoryName = CheckPath('metropolis');
+    MhDirectoryName = CheckPath('metropolis');
 elseif strcmpi(type,'gsa')
-  MhDirectoryName = CheckPath('GSA');
+    MhDirectoryName = CheckPath('GSA');
 else
-  MhDirectoryName = CheckPath('prior');
+    MhDirectoryName = CheckPath('prior');
 end
 if strcmpi(type,'posterior')
-  load([ MhDirectoryName filesep  M_.fname '_mh_history.mat'])
-  TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-  NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+    load([ MhDirectoryName filesep  M_.fname '_mh_history.mat'])
+    TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+    NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
 elseif strcmpi(type,'gsa')
-  load([ MhDirectoryName filesep  M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
-  x=[lpmat0(istable,:) lpmat(istable,:)];
-  clear lpmat istable
-  NumberOfDraws=size(x,1);
-  B=NumberOfDraws; options_.B = B;
+    load([ MhDirectoryName filesep  M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
+    x=[lpmat0(istable,:) lpmat(istable,:)];
+    clear lpmat istable
+    NumberOfDraws=size(x,1);
+    B=NumberOfDraws; options_.B = B;
 else% type = 'prior'
-  NumberOfDraws = 500;
+    NumberOfDraws = 500;
 end
 if ~strcmpi(type,'gsa')
   B = min([round(.5*NumberOfDraws),500]); options_.B = B;
 end
-try delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat'])
-catch disp('No _IRFs (dsge) files to be deleted!')
+try 
+    delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat'])
+catch 
+    disp('No _IRFs (dsge) files to be deleted!')
 end
-try delete([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat'])
-catch disp('No _IRFs (bvar-dsge) files to be deleted!')
+try 
+    delete([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat'])
+catch 
+    disp('No _IRFs (bvar-dsge) files to be deleted!')
 end
 irun = 0;
 IRUN = 0;
@@ -113,14 +118,14 @@ NumberOfIRFfiles_dsgevar = 1;
 ifil2 = 1;
 % Create arrays
 if B <= MAX_nruns
-  stock_param = zeros(B, npar);
+    stock_param = zeros(B, npar);
 else
-  stock_param = zeros(MAX_nruns, npar);
+    stock_param = zeros(MAX_nruns, npar);
 end
 if B >= MAX_nirfs_dsge
-  stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
+    stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
 else
-  stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B);
+    stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B);
 end
 if MAX_nirfs_dsgevar
     if B >= MAX_nirfs_dsgevar
@@ -130,9 +135,9 @@ if MAX_nirfs_dsgevar
     end
     [mYY,mXY,mYX,mXX,Ydata,Xdata] = ...
         var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,...
-                           options_.varlag,-1,options_.datafile,options_.varobs,...
+                           options_.dsge_varlag,-1,options_.datafile,options_.varobs,...
                            options_.xls_sheet,options_.xls_range);
-    NumberOfLags = options_.varlag;
+    NumberOfLags = options_.dsge_varlag;
     NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
     if options_.noconstant
         NumberOfParametersPerEquation = NumberOfLagsTimesNvobs;
@@ -158,30 +163,35 @@ localVars.irun2=irun2;
 localVars.nosaddle=nosaddle;
 
 localVars.type=type;
-if strcmpi(type,'posterior'),
-while b<B 
-    b = b + 1;
-    x(b,:) = GetOneDraw(type);
-end
+if strcmpi(type,'posterior')
+    while b<B 
+        b = b + 1;
+        x(b,:) = GetOneDraw(type);
+    end
 end
 
 if ~strcmpi(type,'prior'),
-   localVars.x=x;
+    localVars.x=x;
 end
 
 b=0;
-
+if options_.dsge_var
+    localVars.gend = gend;
+    localVars.nvobs = nvobs;
+    localVars.NumberOfParametersPerEquation = NumberOfParametersPerEquation;
+    localVars.NumberOfLags = options_.dsge_varlag;
+    localVars.NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
+    localVars.Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
+end 
 localVars.nvar=nvar;
 localVars.IndxVariables=IndxVariables;
 localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
 localVars.MAX_nirfs_dsge=MAX_nirfs_dsge;
 localVars.MAX_nruns=MAX_nruns;
-
 localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
 localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
 localVars.ifil2=ifil2;
 
-
 % Like sequential execution!
 if isnumeric(options_.parallel),% | isunix, % For the moment exclude unix platform from parallel implementation.
     [fout] = PosteriorIRF_core1(localVars,1,B,0);
@@ -349,10 +359,11 @@ localVars.tit=tit;
 localVars.nn=nn;
 localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
 localVars.HPDIRF=HPDIRF;
+localVars.HPDIRFdsgevar=HPDIRFdsgevar;
+localVars.MeanIRFdsgevar = MeanIRFdsgevar;
 localVars.varlist=varlist;
 localVars.MaxNumberOfPlotPerFigure=MaxNumberOfPlotPerFigure;
 
-
 %%% The files .TeX are genereted in sequential way always!
 
 if options_.TeX
@@ -414,8 +425,4 @@ end
 % END parallel code!
 
 
-fprintf('MH: Posterior IRFs, done!\n');
-
-
-
-
+fprintf('MH: Posterior IRFs, done!\n');
\ No newline at end of file

matlab/PosteriorIRF_core1.m

@@ -42,8 +42,6 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,npar,whoiam, ThisMatlab)
 
 global options_ estim_params_ oo_ M_ bayestopt_
 
-
-
 if nargin<4,
     whoiam=0;
 end
@@ -61,6 +59,15 @@ if ~strcmpi(type,'prior'),
     x=myinputs.x;
 end
 
+if options_.dsge_var
+    gend=myinputs.gend;
+    nvobs=myinputs.nvobs;
+    NumberOfParametersPerEquation = myinputs.NumberOfParametersPerEquation;
+    NumberOfLags = myinputs.NumberOfLags;
+    NumberOfLagsTimesNvobs = myinputs.NumberOfLagsTimesNvobs;
+    Companion_matrix = myinputs.Companion_matrix;
+end
+
 nvar=myinputs.nvar;
 IndxVariables=myinputs.IndxVariables;
 MAX_nirfs_dsgevar=myinputs.MAX_nirfs_dsgevar;
@@ -123,23 +130,18 @@ if whoiam
    NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar(whoiam);
 end
 
- while fpar<npar % Parallel 'while'!!!
-   
+while fpar<npar % Parallel 'while'!!!
     fpar = fpar + 1;
     irun = irun+1;
     irun2 = irun2+1;
-    if strcmpi(type,'prior')
-       
-         deep = GetOneDraw(type);
-       
+    if strcmpi(type,'prior')   
+        deep = GetOneDraw(type);
     else
         deep = x(fpar,:);
     end
     stock_param(irun2,:) = deep;  
     set_parameters(deep);
-    [dr,info] = resol(oo_.steady_state,0);
- 
-    
+    [dr,info] = resol(oo_.steady_state,0); 
     if info(1)
         nosaddle = nosaddle + 1;
         fpar = fpar - 1;
@@ -176,8 +178,7 @@ end
     end
     if MAX_nirfs_dsgevar
         IRUN = IRUN+1;
-        %tmp_dsgevar = zeros(options_.irf,nvobs*M_.exo_nbr);
-        [fval,cost_flag,info,PHI,SIGMAu,iXX] =  DsgeVarLikelihood(deep(fpar,:)',gend);
+        [fval,cost_flag,info,PHI,SIGMAu,iXX] =  DsgeVarLikelihood(deep',gend);
         dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
         DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
         SIGMA_inv_upper_chol = chol(inv(SIGMAu*gend*(dsge_prior_weight+1))); 
@@ -194,15 +195,7 @@ end
             explosive_var = any(abs(eig(Companion_matrix))>1.000000001);
         end
         % Get the mean 
-% $$$         if options_.noconstant
-            mu = zeros(1,nvobs); 
-% $$$         else
-% $$$             AA = eye(nvobs);
-% $$$             for lag=1:NumberOfLags
-% $$$                 AA = AA-PHI_draw((lag-1)*nvobs+1:lag*nvobs,:);
-% $$$             end
-% $$$             mu = transpose(AA\transpose(PHI_draw(end,:)));
-% $$$         end
+        mu = zeros(1,nvobs); 
         % Get rotation
         if dsge_prior_weight > 0
             Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
@@ -234,7 +227,7 @@ end
                      int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];,
             eval(['save ' instr]);
             if RemoteFlag==1,
-            OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
+                OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
             end
             NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1; 
             IRUN =0;
@@ -251,14 +244,14 @@ end
                 eval(['save ' instr]);
                 NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
                 if RemoteFlag==1,
-                OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
+                    OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
                 end
                 irun = 0;
             end
         end
         save([MhDirectoryName '/' M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat'],'stock_irf_dsge');
         if RemoteFlag==1,
-        OutputFileName_dsge = [OutputFileName_dsge; {[MhDirectoryName filesep], [M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat']}];
+            OutputFileName_dsge = [OutputFileName_dsge; {[MhDirectoryName filesep], [M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat']}];
         end        
         NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1;
         irun = 0;
@@ -270,9 +263,8 @@ end
         stock = stock_param;
         save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2) '.mat'],'stock');
         if RemoteFlag==1,
-        OutputFileName_param = [OutputFileName_param; {[MhDirectoryName filesep], [M_.fname '_param_irf' int2str(ifil2) '.mat']}];
+            OutputFileName_param = [OutputFileName_param; {[MhDirectoryName filesep], [M_.fname '_param_irf' int2str(ifil2) '.mat']}];
         end
-        
         ifil2 = ifil2 + 1;
         irun2 = 0;
     end
@@ -283,12 +275,11 @@ end
     end
    % if mod(fpar,10)==0 & whoiam,
    if whoiam,
-        fprintf('Done! \n');
-        waitbarString = [ 'Subdraw ' int2str(fpar) '/' int2str(npar) ' done.'];
-        fMessageStatus((fpar-fpar0)/(npar-fpar0),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo)
-    end
-  
- end
+       fprintf('Done! \n');
+       waitbarString = [ 'Subdraw ' int2str(fpar) '/' int2str(npar) ' done.'];
+       fMessageStatus((fpar-fpar0)/(npar-fpar0),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo)
+   end
+end
 
 if whoiam==0
     if nosaddle
@@ -303,10 +294,6 @@ if whoiam==0
     end
 end
 
-
-
-
-
 % Copy the rusults of computation on the call machine (specifically in the
 % directory on call machine that contain the model).
 

matlab/PosteriorIRF_core2.m

@@ -54,6 +54,9 @@ tit=myinputs.tit;
 nn=myinputs.nn;
 MAX_nirfs_dsgevar=myinputs.MAX_nirfs_dsgevar;
 HPDIRF=myinputs.HPDIRF;
+HPDIRFdsgevar=myinputs.HPDIRFdsgevar;
+MeanIRFdsgevar=myinputs.MeanIRFdsgevar;
+
 varlist=myinputs.varlist;
 MaxNumberOfPlotPerFigure=myinputs.MaxNumberOfPlotPerFigure;