From 21f989819e1cb9d5a36d699947896f1523f88c54 Mon Sep 17 00:00:00 2001
From: michel <michel@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Sat, 3 Jun 2006 19:45:05 +0000
Subject: [PATCH] v4: reduced state space for likelihood computation

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@781 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 matlab/DiffuseLikelihood1.m  |  2 +-
 matlab/DiffuseLikelihoodH1.m |  2 +-
 matlab/DsgeLikelihood.m      |  4 ++++
 matlab/DsgeSmoother.m        |  1 +
 matlab/dynare_estimation.m   | 14 +++++++++++++-
 matlab/set_state_space.m     |  8 ++++----
 6 files changed, 24 insertions(+), 7 deletions(-)

diff --git a/matlab/DiffuseLikelihood1.m b/matlab/DiffuseLikelihood1.m
index c20ae2542..8be81e596 100644
--- a/matlab/DiffuseLikelihood1.m
+++ b/matlab/DiffuseLikelihood1.m
@@ -70,7 +70,7 @@ function [LIK, lik] = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,Y,trend,start)
       lik(t)    = log(dF)+transpose(v)*iF*v;
       K         = Pstar(:,mf)*iF; %% premultiplication by the transition matrix T is removed (stephane)
       a         = T*(a+K*v);		%% --> factorization of the transition matrix...
-      Pstar     = T*(Pstar-Pstar(:,mf)*iF*Pstar(mf,:))*transpose(T)+QQ;	%% ... idem (stephane)
+      Pstar     = T*(Pstar-K*Pstar(mf,:))*transpose(T)+QQ;	%% ... idem (stephane)
     end
     notsteady = ~(max(max(abs(Pstar-oldPstar)))<crit);
   end
diff --git a/matlab/DiffuseLikelihoodH1.m b/matlab/DiffuseLikelihoodH1.m
index 41fbc2052..ef6801d46 100644
--- a/matlab/DiffuseLikelihoodH1.m
+++ b/matlab/DiffuseLikelihoodH1.m
@@ -141,7 +141,7 @@ function [LIK, lik] = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,Y,trend,start)
       lik(t)    = log(dF)+transpose(v)*iF*v;
       K         = Pstar(:,mf)*iF; %% premultiplication by the transition matrix T is removed (stephane)
       a         = T*(a+K*v);		%% --> factorization of the transition matrix...
-      Pstar     = T*(Pstar-Pstar(:,mf)*iF*Pstar(mf,:))*transpose(T)+QQ;	%% ... idem (stephane)
+      Pstar     = T*(Pstar-K*Pstar(mf,:))*transpose(T)+QQ;	%% ... idem (stephane)
     end
     notsteady = ~(max(max(abs(Pstar-oldPstar)))<crit);
   end
diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m
index 857857b0d..a6bad363a 100644
--- a/matlab/DsgeLikelihood.m
+++ b/matlab/DsgeLikelihood.m
@@ -87,6 +87,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
   % 2. call model setup & reduction program
   %------------------------------------------------------------------------------
   [T,R,SteadyState,info] = dynare_resolve;
+  rs = bayestopt_.restrict_state;
   if info(1) == 1 | info(1) == 2 | info(1) == 5
     fval = bayestopt_.penalty+1;
     cost_flag = 0;
@@ -96,6 +97,9 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
     cost_flag = 0;
     return
   end
+  T = T(rs,rs);
+  R = R(rs,:);
+  bayestopt_.mf = bayestopt_.mf1;
   if options_.loglinear == 1
     constant = log(SteadyState(bayestopt_.mfys));
   else
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 3e40f1110..2a6655294 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -16,6 +16,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoo
   % 2. call model setup & reduction program
   %------------------------------------------------------------------------------
   [T,R,SteadyState] = dynare_resolve;
+  bayestopt_.mf = bayestopt_.mf2;
   if options_.loglinear == 1
     constant = log(SteadyState(bayestopt_.mfys));
   else
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index a23e248e4..e90b91f65 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -200,8 +200,20 @@ for i=1:n_varobs
   k = [k strmatch(deblank(options_.varobs(i,:)),M_.endo_names(dr.order_var,:),'exact')];
   k1 = [k1 strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')];
 end
+% union of observed and state variables
+k2 = union(k',[dr.nstatic+1:dr.nstatic+dr.npred]');
+% including variables in t-2 and earlier, if any
+k2 = [k2;[M_.endo_nbr+(1:dr.nspred-dr.npred)]'];
 
-bayestopt_.mf 	= k;
+% set restrict_state to postion of observed + state variables
+% in expanded state vector
+bayestopt_.restrict_state = k2;
+% set mf1 to positions of observed variables in restricted state vector
+% for likelihood computation
+[junk,bayestopt_.mf1] = ismember(k,k2); 
+% set mf2 to positions of observed variables in expanded state vector
+% for filtering and smoothing
+bayestopt_.mf2 	= k;
 bayestopt_.mfys = k1;
 options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
 gend = options_.nobs;
diff --git a/matlab/set_state_space.m b/matlab/set_state_space.m
index 6ed41cced..5293a162b 100644
--- a/matlab/set_state_space.m
+++ b/matlab/set_state_space.m
@@ -78,7 +78,7 @@ dr.kad = kad;
 dr.kae = kae;
 dr.nboth = nboth;
 dr.nfwrd = nfwrd;
-
-% 06/01/03 MJ removed negative entries for kstate(:,3) when lags on more
-%             than 1 period
-% 06/25/03 MJ corrected logic for kmask computation when M_.maximum_lead == 0
+% number of forward variables in the state vector
+dr.nsfwrd = sum(kstate(:,2) > M_.maximum_lag+1);
+% number of predetermined variables in the state vector
+dr.nspred = sum(kstate(:,2) <= M_.maximum_lag+1);