Corrected bug related to the constant in bvar-dsge.

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1484 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/DsgeVarLikelihood.m

@@ -34,6 +34,9 @@ ns = nvx+nvn+ncx+ncn;
 NumberOfObservedVariables = size(options_.varobs,1);
 NumberOfLags = options_.varlag;
 NumberOfParameters = NumberOfObservedVariables*NumberOfLags ;
+if ~options_.noconstant
+    NumberOfParameters = NumberOfParameters + 1;
+end
 
 mYY = evalin('base', 'mYY');
 mYX = evalin('base', 'mYX');
@@ -108,8 +111,7 @@ if ~options_.noconstant
         constant = transpose(log(SteadyState(bayestopt_.mfys)));
     else
         constant = transpose(SteadyState(bayestopt_.mfys));
-    end
-    NumberOfParameters = NumberOfParameters + 1;	
+    end	
 else
     constant = zeros(1,NumberOfObservedVariables);
 end
@@ -153,7 +155,7 @@ for i = 1:NumberOfLags-1
 end
 if ~options_.noconstant
     % Add one row and one column to GXX
-    GXX = [GXX , ones(NumberOfLags*NumberOfObservedVariables,1) ; ones(1,NumberOfParameters)];
+    GXX = [GXX , kron(ones(NumberOfLags,1),constant') ; [  kron(ones(1,NumberOfLags),constant) , 1] ];
 end
 
 GYY = TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1);

matlab/PosteriorIRF.m

@@ -123,6 +123,11 @@ if MAX_nirfs_dsgevar
         var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,options_.varlag,-1);
     NumberOfLags = options_.varlag;
     NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
+    if options_.noconstant
+        NumberOfParametersPerEquation = NumberOfLagsTimesNvobs;
+    else
+        NumberOfParametersPerEquation = NumberOfLagsTimesNvobs+1;
+    end
     Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
 end
 b = 0;
@@ -175,7 +180,7 @@ while b<=B
     end
     if MAX_nirfs_dsgevar
         IRUN = IRUN+1;
-        tmp_dsgevar = zeros(options_.irf,nvobs*M_.exo_nbr);
+        %tmp_dsgevar = zeros(options_.irf,nvobs*M_.exo_nbr);
         [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));
@@ -183,15 +188,25 @@ while b<=B
         explosive_var  = 1;
         while explosive_var
             % draw from the marginal posterior of SIGMA
-            SIGMAu_draw = rand_inverse_wishart(nvobs, DSGE_PRIOR_WEIGHT-NumberOfLagsTimesNvobs, ...
+            SIGMAu_draw = rand_inverse_wishart(nvobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ...
                                                SIGMA_inv_upper_chol);
             % draw from the conditional posterior of PHI
-            PHI_draw = rand_matrix_normal(NumberOfLagsTimesNvobs,nvobs, PHI, ...
+            PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,nvobs, PHI, ...
                                            chol(SIGMAu_draw)', chol(iXX)');
-            Companion_matrix(1:nvobs,:) = transpose(PHI_draw);
+            Companion_matrix(1:nvobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:));
             % Check for stationarity
             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
         % Get rotation
         if dsge_prior_weight > 0
             Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
@@ -207,8 +222,9 @@ while b<=B
         for t = 2:options_.irf
             PHIpower = Companion_matrix*PHIpower;
             tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
-            irfs(t,:)  = tmp3(:)';
-        end            
+            irfs(t,:)  = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu);
+        end
+        tmp_dsgevar = kron(ones(options_.irf,1),mu);
         for j = 1:(nvobs*M_.exo_nbr)
             if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10 
                 tmp_dsgevar(:,j) = (irfs(:,j));