From c14493361b7190ed9082296df0276a4e0630ab41 Mon Sep 17 00:00:00 2001
From: adjemian <adjemian@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Wed, 14 Nov 2007 17:42:18 +0000
Subject: [PATCH] Added estimation of bvar-dsge models with constant.

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1441 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 matlab/DsgeVarLikelihood.m  | 48 +++++++++++++++++++++++--------------
 matlab/dynare_estimation.m  |  4 ----
 matlab/var_sample_moments.m |  8 +++----
 3 files changed, 34 insertions(+), 26 deletions(-)

diff --git a/matlab/DsgeVarLikelihood.m b/matlab/DsgeVarLikelihood.m
index 2dc5f0cf9..a096f2984 100644
--- a/matlab/DsgeVarLikelihood.m
+++ b/matlab/DsgeVarLikelihood.m
@@ -14,8 +14,7 @@ function [fval,cost_flag,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(xparam1,gend)
 %   o iXX           [double]     inv(X'X).
 %
 % ALGORITHM
-%   None.       
-%
+%   None.
 % SPECIAL REQUIREMENTS
 %   None.
 %  
@@ -83,7 +82,7 @@ dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
 % Is the DSGE prior proper?
 if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/gend;
     fval = bayestopt_.penalty*min(1e3,(NumberOfParameters+NumberOfObservedVariables)/gend-dsge_prior_weight);
-    info = 51
+    info = 51;
     cost_flag = 0;
     return;
 end
@@ -104,10 +103,15 @@ elseif info(1) == 3 | info(1) == 4 | info(1) == 20
     cost_flag = 0;
     return
 end
-if options_.loglinear == 1
-    constant = log(SteadyState(bayestopt_.mfys));
+if ~options_.noconstant
+    if options_.loglinear 
+        constant = transpose(log(SteadyState(bayestopt_.mfys)));
+    else
+        constant = transpose(SteadyState(bayestopt_.mfys));
+    end
+    NumberOfParameters = NumberOfParameters + 1;	
 else
-    constant = SteadyState(bayestopt_.mfys);
+    constant = zeros(1,NumberOfObservedVariables);
 end
 if bayestopt_.with_trend == 1
     disp('DsgeVarLikelihood :: Linear trend is not yet implemented!')
@@ -115,33 +119,41 @@ if bayestopt_.with_trend == 1
 end
 
 %------------------------------------------------------------------------------
-% 3. theorretical moments (second order)
+% 3. theoretical moments (second order)
 %------------------------------------------------------------------------------
 tmp0 = lyapunov_symm(T,R*Q*R');% I compute the variance-covariance matrix
-                               % of the restricted state vector.
-%bayestopt_.mf = bayestopt_.mf1;????????
-mf  = bayestopt_.mf1;
-
+mf  = bayestopt_.mf1;          % of the restricted state vector.
 
+% Get the non centered second order moments
 TheoreticalAutoCovarianceOfTheObservedVariables = ...
     zeros(NumberOfObservedVariables,NumberOfObservedVariables,NumberOfLags+1);
-TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) = tmp0(mf,mf);
+TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) = tmp0(mf,mf)+constant'*constant;
 for lag = 1:NumberOfLags
   tmp0 = T*tmp0;
-  TheoreticalAutoCovarianceOfTheObservedVariables(:,:,lag+1) = tmp0(mf,mf);
+  TheoreticalAutoCovarianceOfTheObservedVariables(:,:,lag+1) = tmp0(mf,mf) ...
+      + constant'*constant;
 end
+% Build the theoretical "covariance" between Y and X
 GYX = zeros(NumberOfObservedVariables,NumberOfParameters);
 for i=1:NumberOfLags
   GYX(:,(i-1)*NumberOfObservedVariables+1:i*NumberOfObservedVariables) = ...
       TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1);
 end
+if ~options_.noconstant
+    GYX(:,end) = constant';
+end
+% Build the theoretical "covariance" between X and X
 GXX = kron(eye(NumberOfLags), ...
            TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1));
 for i = 1:NumberOfLags-1
-  tmp1 = diag(ones(NumberOfLags-i,1),i); 
-  tmp2 = diag(ones(NumberOfLags-i,1),-i);
-  GXX = GXX + kron(tmp1,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1));
-  GXX = GXX + kron(tmp2,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1)');
+    tmp1 = diag(ones(NumberOfLags-i,1),i); 
+    tmp2 = diag(ones(NumberOfLags-i,1),-i);
+    GXX = GXX + kron(tmp1,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1));
+    GXX = GXX + kron(tmp2,TheoreticalAutoCovarianceOfTheObservedVariables(:,:,i+1)');
+end
+if ~options_.noconstant
+    % Add one row and one column to GXX
+    GXX = [GXX , ones(NumberOfLags*NumberOfObservedVariables,1) ; ones(1,NumberOfParameters)];
 end
 
 GYY = TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1);
@@ -150,7 +162,7 @@ assignin('base','GYY',GYY);
 assignin('base','GXX',GXX);
 assignin('base','GYX',GYX);
 
-if ~isinf(dsge_prior_weight) 
+if ~isinf(dsge_prior_weight)
   tmp0 = dsge_prior_weight*gend*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ;
   tmp1 = dsge_prior_weight*gend*GYX + mYX;
   tmp2 = inv(dsge_prior_weight*gend*GXX+mXX);
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index b581dd1a7..1ef97779d 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -25,9 +25,6 @@ end
 if options_.prefilter == 1
   options_.noconstant = 1;
 end
-if options_.bvar_dsge 
-  options_.noconstant = 1;
-end
 
 if options_.filtered_vars ~= 0 & options_.filter_step_ahead == 0
   options_.filter_step_ahead = 1;
@@ -223,7 +220,6 @@ end
 
 %% compute sample moments if needed (bvar-dsge)
 if options_.bvar_dsge~isempty(strmatch('dsge_prior_weight',M_.param_names))
-    options_.noconstant = 1;
     if options_.noconstant
         evalin('base',['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ...
                  'var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,options_.varlag,-1);'])
diff --git a/matlab/var_sample_moments.m b/matlab/var_sample_moments.m
index 2b6fa7704..f4a2bea72 100644
--- a/matlab/var_sample_moments.m
+++ b/matlab/var_sample_moments.m
@@ -82,7 +82,7 @@ if qlag > FirstObservation
 end
 
 NumberOfObservations = LastObservation-FirstObservation+1;% This is T.
-NumberOfVariables = options_.varobs;% This is m.
+NumberOfVariables = size(options_.varobs,1);% This is m.
 if var_trend_order == -1% No constant no linear trend case.
     X = zeros(NumberOfObservations,NumberOfVariables*qlag);
 elseif var_trend_order == 0% Constant and no linear trend case.
@@ -101,12 +101,12 @@ Y = data(FirstObservation:LastObservation,:);
 for t=1:NumberOfObservations
   line = t + FirstObservation-1;
   for lag = 1:qlag
-    X(t,(lag-1)*NumberOfVariables+1:lag*NumberOfVariables) = data(line-lag,:);
+      X(t,(lag-1)*NumberOfVariables+1:lag*NumberOfVariables) = data(line-lag,:);
   end
   if var_trend_order == 0
-    X(t,indx) = ones(1,NumberOfVariables);
+      X(t,indx) = 1;
   elseif var_trend_order == 1
-  X(t,indx) = [ 1 , t ];
+      X(t,indx) = [ 1 , t ];
   end
 end