function [StateVarStructure,status,kalmanstatus] = CalculateDSGEStateVariables(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,ObsVarDec)
% CalculateDSGEStateVariables: Calculates the forecasted, updated and smoothed state variables as well as
%                              the forecasted observed variables
%
% USAGE:
%
%       [StateVarStructure,status,kalmanstatus] = CalculateDSGEStateVariables(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,ObsVarDec)
%
% REQUIRED INPUT:  theta (vector) with values of the original DSGE model parameters
%
%                  thetaPositions (structure) with fields parameter and distribution. This function
%                                 needs the former field to know which parameter name corresponds to
%                                 position i in theta.
%
%                  ModelParameters (structure) whose fields have names given by the parameter names in the
%                                 DSGE model.
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure) with initialization information.
%
%                  ObsVarDec (boolean) which determines if the observed variables decomposition is calculated
%                                 or not.
%
% REQUIRED OUTPUT: StateVarStructure (structure) with fields "Y", "Yhat" containing the observed variables and
%                                 the 1-step ahead forecasts of these, "Ksitt1", "Ksitt", and "KsitT" with
%                                 the forecasted, updated, and smoothed values for the state variables. Furthermore,
%                                 the field "MaxEigenvalue" contains the largest eigenvalue of the state transition
%                                 matrix F, "etatT" the smooth projections of the unique state shocks, "wtT" the
%                                 smooth projections of the measurement errors, and "wtNames" the names of the
%                                 measurement errors (determined directly from the observed variable's names. The
%                                 field "wthT" contains all the measurement errors. Furthermore, "TrainSample" is
%                                 a boolean variable that indicates if a training sample was employed for the
%                                 Kalman filter. The field "X" holds the exogneous variables, while the field
%                                 "XiInit" contains the matrix with typical element H'*(F^t)*Ksi0T that is used
%                                 for observed variable decompositions. Similarly, the field "etaDecomp" contains
%                                 a 3D matrix with the parts of the observed variable decomposition that are explained
%                                 by the q economic shocks. The field "KeepVar" contains the column number of B0
%                                 that are non-zero. The parameter matrices of the measurement equation are given by
%                                 the fields "A", "H", and "R", while the fields "F" a "B0" contain the parameter
%                                 matrices of the state equation.
%
% OPTIONAL OUTPUT: status (integer) Indicates whether the solution of the model is unique or not.
%                                 mcode is 1 when unique, 0 if not computed, and various other values when
%                                 there are problems with the solution. Specifically:
%
%                                      mcode = 2:      Roots not correctly computed by real_schur
%                                      mcode = 3:      Too many big roots
%                                      mcode = 35:     Too many big roots, and q(:,right) is singular
%                                      mcode = 4:      Too few big roots
%                                      mcode = 45:     Too few big roots, and q(:,right) is singular
%                                      mcode = 5:      q(:,right) is singular
%                                      mcode = 61:     Too many exact shiftrights
%                                      mcode = 62:     Too many numeric shiftrights
%                                      mcode = 7:      The a matrix has infinite or NaN values. Eigenvalues
%                                                      cannot be calculated.
%                                      mcode = 8:      The function "compute_aim_matrices" returns complex
%                                                      numbers.
%                                      else            Return code not properly specified
%
%                  kalmanstatus (boolean) indicating the covariance stationarity of the state variables. If all the
%                                 eigenvalues of F are inside the unit circle, then kalmanstatus=0 and the log-likelihood function
%                                 is calculated. If not, then status=1 and initP is set to 0.
%
% NOTE: Three output arguments should only be used when theta has the initial values.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2009 European Central Bank.
%
%                       First version: August 21, 2006.
%                        This version. March 9, 2009.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 22-08-2006: Added the try, catch code around the StateSmoother function in case the reduced rank fixes don't work. the
%               option DSGEModel.KalmanFirstObservation was also added among the KalmanFilter function inputs.
%
% * 25-08-2006: The compute_aim_matrices, measurement equation file, and the file with parameters to update are now
%               always copied to the tmp directory to make sure that the correct version is located there.
%
% * 28-08-2006: Made sure that the StateVarStructure.Y output takes the variable DSGEModel.KalmanFirstObservation into account.
%
% * 14-09-2006: Added the changes to the measurement equation call.
%
% * 11-10-2006: Added computations of state shocks.
%
% * 12-10-2006: Added computation of the measurement errors.
%
% * 16-10-2006: Fixed the output for the measurement errors such that only the measurement errors with non-zero variance
%               and sent as output. The equation information about the measurement errors is now also given.
%
% * 18-10-2006: Made sure that the B0 matrix used when estimating the state shocks has full column rank.
%
% * 31-10-2006: Added the input "DAToleranceValue" for the call to the KalmanFilter function.
%
% * 10-11-2006: Added A, H, F, and B0 to the StateVarStructure output.
%
% * 14-11-2006: Included the field KeepVar in the output structure. Moreover, the optional output P01 from KalmanFilter is
%               takes into account for possibly calculation the smoothed values of the state vector at time 0. This requires
%               that a training sample is used in the Kalman filter calculations.
%
% * 15-11-2006: Fixed a bug in the call to the StateSmoother function where the training sample was wrongly included for
%               the observed variables.
%
% * 16-11-2006: Added the field "R" to the StateVarStructure structure. The value gives the covariance matrix of the
%               measurement errors.
%
% * 17-11-2006: Added the field "lnLt" to the StateVarStructure structure.
%
% * 20-11-2006: The input vector InitialStateVector for the Kalman filter is now created based on user preferences
%               and sent to the function.
%
% * 23-11-2006: Added try-catch statements around eval calls.
%
% * 04-12-2006: Added the "XiInit" and "etaDecomp" fields for the observed variable decomposition.
%
% * 05-12-2006: Updated the documentation about the StateVarStructure fields.
%
% * 06-12-2006: Added the input variable "ObsVarDec".
%
% * 04-01-2007: Made changes such that the full sample is used for all computations, while the estimation sample
%               is used for setting up the output. The output structure now includes "XietaDecomp" and "XietaInit"
%               with historical decompositions of the smoothed estimates of the state variables, and initialization
%               values, respectively.
%
% * 20-02-2007: The function now also tries to estimate the updated economic shocks and measurement errors.
%
% * 09-03-2007: Made sure that the missing input "CurrINI" was included.
%
% * 27-03-2007: Updated the documentation and made some output checks.
%
% * 13-11-2007: Updated the documentation.
%
% * 03-12-2007: Made sure the Kalman filter function includes the correct X matrix.
%
% * 21-12-2007: Added the code for dealing with unit-root state variables when running the Kalman filter.
%
% * 16-04-2008: The code now checks if the model has a unique convergent solution also at the posterior mode
%               estimates.
%
% * 07-05-2008: Added error catching code for the Kalman filter and error display for the filter and the
%               smoother.
%
% * 16-05-2008: Took the "DSGEModel.AIMTolerance" input for AiMSolver into account.
%
% * 23-05-2008: Updated the documentation.
%
% * 27-05-2008: The mcode value 7 is now supported.
%
% * 02-07-2008: Added the variable "DSGEModel.AllowUnitRoot" for the Kalman filters.
%
% * 03-07-2008: Made sure that the measurement equation receives only the active state variable names.
%               Added the try-catch code for AiMSolver.
%
% * 07-07-2008: Added code to allow for the case when the measurement matrix H is time-varying.
%
% * 09-07-2008: Made sure that a time-varying measurement matrix H is exported for the correct
%               sample.
%
% * 22-01-2009: Updated the documentation.
%
% * 30-01-2009: Made sure that status is only set to mcode when the latter exists.
%
% * 09-02-2009: Changed the input "length(DSGEModel.StateVariablePositions)" to "size(F,1)" when calling the
%               measurement equation function.
%
% * 13-02-2009: The mcode values 7 and 8 are now properly mentioned when encountered.
%
% * 09-03-2009: Added quotation marks around paths for the dos copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize the output
%
StateVarStructure = [];
if nargin<6;
   ObsVarDec = 0;
end;
%
% set up the ModelParameters structure based on theta and thetaPositions
%
if nargout==1;
   %
   % Make sure that the ModelParameter structure takes the correct values.
   % This is not necessary to do when initial values are given to the function.
   % In that case, more substantial checking is required!
   %
   ModelParameters = ThetaToModelParameters(ModelParameters,theta,thetaPositions);
   %
   % copy files to the tmp folder
   %
   [stat,msg] = dos(['copy /Y "' GetPath(DSGEModel.AIMDataFile) 'compute_aim_matrices.m" "' pwd '\tmp"']);
   [stat,msg] = dos(['copy /Y "' DSGEModel.MeasurementEquationFile '" "' pwd '\tmp"']);
   if FileExist(DSGEModel.UpdateParameterFile)==1;
      [stat,msg] = dos(['copy /Y "' DSGEModel.UpdateParameterFile '" "' pwd '\tmp"']);
   end;
   %
   % Update any parameteters that need updating
   %
   if FileExist(DSGEModel.UpdateParameterFile)==1;
      eval(['ModelParameters = ' GetName(DSGEModel.UpdateParameterFile) '(ModelParameters);']);
   end;
end;
%
% load the AIMData file
%
ErrorStr = '';
try;
   AIMData = load(DSGEModel.AIMDataFile);
catch;
   ErrorStr = ['Unable to load the file "' DSGEModel.AIMDataFile '". Message caught is: ' lasterr];
end;
if isempty(ErrorStr)==0;
   if nargout>1;
      status = 0;
      kalmanstatus = 0;
   end;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% See if we can solve the model
%
if nargout==1;
   [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
else;
   try;
      [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
   catch;
      if exist('mcode')>0;
         status = mcode;
      else;
         status = 0;
      end;
      ErrorStr = ['There was an error running AiM. YADA received the following message: ' lasterr];
      About(ErrorStr,'error','AiM Error',300,500,CurrINI);
      drawnow;
      return;
   end;
end;
if nargout>2;
   status = mcode;
   if status~=1;
      return;
   end;
else;
   if mcode~=1;
      if mcode==2;
         mcodeStr = 'Roots are not correctly computed by real_schur.'
      elseif mcode==3;
         mcodeStr = 'Too many big roots.';
      elseif mcode==35;
         mcodeStr = 'Too many big roots, and q(:,right) is singular.';
      elseif mcode==4;
         mcodeStr = 'Too few big roots.';
      elseif mcode==45;
         mcodeStr = 'Too few big roots, and q(:,right) is singular.';
      elseif mcode==5;
         mcodeStr = 'q(:,right) is singular.';
      elseif mcode==61;
         mcodeStr = 'Too many exact shiftrights.';
      elseif mcode==62;
         mcodeStr = 'Too many numeric shiftrights.';
      elseif mcode==7;
         mcodeStr = 'Infinite or NaN values detected.';
      elseif mcode==8;
         mcodeStr = 'The function "compute_aim_matrices" returns complex numbers.';
      else;
         mcodeStr = 'Return code not properly specified.';
      end;
      txt = ['The AiM solver provided the return code: ' num2str(mcode) ', i.e., "' mcodeStr '"'];
      About(txt,'information','YADA - AiM Solver Error',200,500,CurrINI);
      return;
   end;
end;
%
% Convert to the state equation matrices
%
[F,B0] = AiMtoStateSpace(BMatrix,SMatrix,DSGEModel.StateVariablePositions,DSGEModel.StateShockPositions);
%
% obtain the measurement equation matrices
%
try;
   eval(['[A,H,R] = ' GetName(DSGEModel.MeasurementEquationFile) '(ModelParameters,AIMData.endog(DSGEModel.StateVariablePositions,:),DSGEModel.VariableNames,DSGEModel.XVariableNames,DSGEModel.n,size(F,1),DSGEModel.k);']);
catch;
   %
   % display an error message
   %
   ErrorStr = ['Error found in "' GetName(DSGEModel.MeasurementEquationFile) '.m". Message caught is: ' lasterr];
   About(ErrorStr,'error','Error - Measurement Equation',150,500,CurrINI);
   return;
end;
if length(size(H))==3;
   if size(H,3)<size(DSGEModel.Y,2);
      %
      % too few data points in the 3D H matrix.
      %
      ErrorStr = ['The number of time periods for the time-varying measurement matrix H (' int2str(size(H,3)) ') is lower than the number of observations (T = ' int2str(size(DSGEModel.Y,2)) '). YADA has therefore aborted from the calculations.'];
      About(ErrorStr,'error','Error - Measurement Equation',200,500,CurrINI);
      return;
   end;
end;
%
% set the InitialStateVector value
%
if DSGEModel.UseOwnInitialState==1;
   if length(DSGEModel.InitialStateValues)==size(AIMData.endog,1);
      DSGEModel.InitialStateVector = DSGEModel.InitialStateValues(DSGEModel.StateVariablePositions);
   else;
      DSGEModel.InitialStateVector = zeros(length(DSGEModel.StateVariablePositions),1);
   end;
else;
   DSGEModel.InitialStateVector = zeros(length(DSGEModel.StateVariablePositions),1);
end;
%
% Get the sub-sample
%
[DSGEModel.FirstPeriod,DSGEModel.LastPeriod] = CreateSubSample(DSGEModel);
%
% We can now try to run the Kalman filter
%
try;
   if length(size(H))==2;
      if isempty(DSGEModel.UnitRootStates)==1;
         if isempty(DSGEModel.X)==0;
            [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
         else;
            [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
         end;
      else;
         if isempty(DSGEModel.X)==0;
            [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
         else;
            [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
         end;
      end;
   else;
      if isempty(DSGEModel.UnitRootStates)==1;
         if isempty(DSGEModel.X)==0;
            [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),A,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
         else;
            [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
         end;
      else;
         if isempty(DSGEModel.X)==0;
            [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),A,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
         else;
            [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = UnitRootKalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,1,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
         end;
      end;
   end;
catch;
   txt = ['YADA caught the following error message from the Kalman filter routine: ' lasterr];
   About(txt,'error','Error - Kalman Filter',200,500,CurrINI);
   return;
end;
%
% we should adjust everything based on DSGEModel.KalmanFirstObservation
%
if nargout>1;
   kalmanstatus = stat;
   if kalmanstatus~=0;
      return;
   end;
end;
%
% Run the state smoother
%
try;
   if length(size(H))==2;
      [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmoother(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),Yhat,Ksitt1,Ptt1,H,F,R);
   else;
      [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmootherHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),Yhat,Ksitt1,Ptt1,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,R);
   end;
catch;
   Ksitt = [];
   KsitT = [];
   txt = ['YADA caught the following error message from the Kalman smoother: ' lasterr];
   About(txt,'error','Error - Kalman Smoother',200,500,CurrINI);
   return;
end;
if DSGEModel.KalmanFirstObservation>1;
   StateVarStructure.TrainSample = 1;
else;
   StateVarStructure.TrainSample = 0;
end;
%
% now we can calculate the shocks
%
[S,KeepVar] = RemoveRedundantColumns(B0);
if isempty(S)==0;
   etatT = inv(S'*S)*S'*(KsitT(:,2:size(KsitT,2))-(F*KsitT(:,1:size(KsitT,2)-1)));
   %
   % as the first observation we assume that xi0T = 0
   %
   etatT = [inv(S'*S)*S'*KsitT(:,1) etatT];
   etatt = inv(S'*S)*S'*(Ksitt(:,2:size(Ksitt,2))-(F*Ksit1t));
   etatt = [inv(S'*S)*S'*Ksitt(:,1) etatt];
   if ObsVarDec==1;
      %
      % need to speed this up
      %
      StateVarStructure.etaDecomp = zeros(size(Yhat,1),size(Yhat,2),size(etatT,1));
      StateVarStructure.XietaDecomp = zeros(size(KsitT,1),size(KsitT,2),size(etatT,1));
      %
      % initialize
      %
      Fi = eye(size(F,1));
      T = length(Yhat);
      for ti=0:T-1;
         %
         % select a shock
         %
         for j=1:size(etatT,1);
            for t=ti+1:T;
               StateVarStructure.XietaDecomp(:,t,j) = StateVarStructure.XietaDecomp(:,t,j)+(Fi*S(:,j)*etatT(j,t-ti));
            end;
         end;
         %
         % update Fi. This is the speediest algorithm since F may have high dimension making
         % F^i a costly computation
         %
         Fi = Fi*F;
         drawnow;
      end;
      if length(size(H))==2;
         for j=1:size(etatT,1);
            StateVarStructure.etaDecomp(:,:,j) = H'*StateVarStructure.XietaDecomp(:,:,j);
         end;
      else;
         for j=1:size(etatT,1);
            StateVarStructure.etaDecomp(:,:,j) = H(:,:,DSGEModel.FirstPeriod+j-1)'*StateVarStructure.XietaDecomp(:,:,j);
         end;
      end;
   else;
      StateVarStructure.XiInit = [];
      StateVarStructure.etaDecomp = [];
      StateVarStructure.XietaDecomp = [];
      StateVarStructure.XietaInit = [];
   end;
else;
   etatT = [];
   etatt = [];
   StateVarStructure.XiInit = [];
   StateVarStructure.etaDecomp = [];
   StateVarStructure.XietaDecomp = [];
   StateVarStructure.XietaInit = [];
end;
%
T = size(Yhat,2);
StateVarStructure.KeepVar = KeepVar;
StateVarStructure.Y = DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod);
StateVarStructure.X = DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod);
StateVarStructure.A = A;
StateVarStructure.H = H;
StateVarStructure.F = F;
StateVarStructure.B0 = B0;
StateVarStructure.R = R;
StateVarStructure.Yhat = Yhat;
StateVarStructure.Ksitt1 = Ksitt1;
StateVarStructure.Ksitt = Ksitt;
StateVarStructure.KsitT = KsitT;
StateVarStructure.etatT = etatT;
StateVarStructure.etatt = etatt;
StateVarStructure.MaxEigenvalue = max(abs(eig(F)));
StateVarStructure.lnLt = lnLt;
%
% Now we estimate the measurement errors based on the full sample information
%
if rank(R)>0;
   if length(size(H))==2;
      if min(size(StateVarStructure.X))>=1;
         wtT = StateVarStructure.Y-(A'*StateVarStructure.X)-(H'*KsitT);
         wtt = StateVarStructure.Y-(A'*StateVarStructure.X)-(H'*Ksitt);
      else;
         wtT = StateVarStructure.Y-(H'*KsitT);
         wtt = StateVarStructure.Y-(H'*Ksitt);
      end;
      StateVarStructure.wtT = [];
      StateVarStructure.wtt = [];
      StateVarStructure.wtNames = '';
      for i=1:length(R);
         if R(i,i)>0;
            StateVarStructure.wtT = [StateVarStructure.wtT;wtT(i,:)];
            StateVarStructure.wtt = [StateVarStructure.wtt;wtt(i,:)];
            StateVarStructure.wtNames = strvcat(StateVarStructure.wtNames,['Error in ' StringTrim(DSGEModel.VariableNames(i,:)) ' equation']);
         end;
      end;
      StateVarStructure.wthT = wtT;
      StateVarStructure.wtht = wtt;
   else;
      wtT = zeros(size(StateVarStructure.Y));
      wtt = wtT;
      for t=1:size(StateVarStructure.Y,2);
         if min(size(StateVarStructure.X))>=1;
            wtT(:,t) = StateVarStructure.Y(:,t)-(A'*StateVarStructure.X(:,t))-(H(:,:,DSGEModel.FirstPeriod+t-1)'*KsitT(:,t));
            wtt(:,t) = StateVarStructure.Y(:,t)-(A'*StateVarStructure.X(:,t))-(H(:,:,DSGEModel.FirstPeriod+t-1)'*Ksitt(:,t));
         else;
            wtT(:,t) = StateVarStructure.Y(:,t)-(H(:,:,DSGEModel.FirstPeriod+t-1)'*KsitT(:,t));
            wtt(:,t) = StateVarStructure.Y(:,t)-(H(:,:,DSGEModel.FirstPeriod+t-1)'*Ksitt(:,t));
         end;
      end;
      StateVarStructure.wtT = [];
      StateVarStructure.wtt = [];
      StateVarStructure.wtNames = '';
      for i=1:length(R);
         if R(i,i)>0;
            StateVarStructure.wtT = [StateVarStructure.wtT;wtT(i,:)];
            StateVarStructure.wtt = [StateVarStructure.wtt;wtt(i,:)];
            StateVarStructure.wtNames = strvcat(StateVarStructure.wtNames,['Error in ' StringTrim(DSGEModel.VariableNames(i,:)) ' equation']);
         end;
      end;
      StateVarStructure.wthT = wtT;
      StateVarStructure.wtht = wtt;
   end;
else;
   StateVarStructure.wtT = [];
   StateVarStructure.wthT = [];
   StateVarStructure.wtt = [];
   StateVarStructure.wtht = [];
   StateVarStructure.wtNames = '';
end;
if ObsVarDec==1;
   %
   % comoute the initial state effect
   %
   etaDecomp = zeros(size(Yhat));
   StateVarStructure.XietaInit = StateVarStructure.KsitT;
   for j=1:size(StateVarStructure.etaDecomp,3);
      etaDecomp = etaDecomp+StateVarStructure.etaDecomp(:,:,j);
      StateVarStructure.XietaInit = StateVarStructure.XietaInit-StateVarStructure.XietaDecomp(:,:,j);
   end;
   etaDecomp = StateVarStructure.Y-(StateVarStructure.A'*StateVarStructure.X)-etaDecomp;
   if isempty(StateVarStructure.wthT)==0;
      etaDecomp = etaDecomp-StateVarStructure.wthT;
   end;
   StateVarStructure.XiInit = etaDecomp;
end;
%
if DSGEModel.KalmanFirstObservation>1;
   %
   % adjust the samples
   %
   StateVarStructure.Y = StateVarStructure.Y(:,DSGEModel.KalmanFirstObservation:T);
   StateVarStructure.X = StateVarStructure.X(:,DSGEModel.KalmanFirstObservation:T);
   StateVarStructure.Yhat = StateVarStructure.Yhat(:,DSGEModel.KalmanFirstObservation:T);
   StateVarStructure.Ksitt1 = StateVarStructure.Ksitt1(:,DSGEModel.KalmanFirstObservation:T);
   StateVarStructure.Ksitt = StateVarStructure.Ksitt(:,DSGEModel.KalmanFirstObservation:T);
   StateVarStructure.KsitT = StateVarStructure.KsitT(:,DSGEModel.KalmanFirstObservation:T);
   StateVarStructure.etatT = StateVarStructure.etatT(:,DSGEModel.KalmanFirstObservation:T);
   StateVarStructure.etatt = StateVarStructure.etatt(:,DSGEModel.KalmanFirstObservation:T);
   StateVarStructure.lnLt = StateVarStructure.lnLt(:,DSGEModel.KalmanFirstObservation:T);
   if isempty(StateVarStructure.wthT)==0;
      StateVarStructure.wthT = StateVarStructure.wthT(:,DSGEModel.KalmanFirstObservation:T);
      StateVarStructure.wtT = StateVarStructure.wtT(:,DSGEModel.KalmanFirstObservation:T);
      StateVarStructure.wtht = StateVarStructure.wtht(:,DSGEModel.KalmanFirstObservation:T);
      StateVarStructure.wtt = StateVarStructure.wtt(:,DSGEModel.KalmanFirstObservation:T);
   end;
   if ObsVarDec==1;
      StateVarStructure.etaDecomp = StateVarStructure.etaDecomp(:,DSGEModel.KalmanFirstObservation:T,:);
      StateVarStructure.XiInit = StateVarStructure.XiInit(:,DSGEModel.KalmanFirstObservation:T);
      StateVarStructure.XietaDecomp = StateVarStructure.XietaDecomp(:,DSGEModel.KalmanFirstObservation:T,:);
      StateVarStructure.XietaInit = StateVarStructure.XietaInit(:,DSGEModel.KalmanFirstObservation:T);
   end;
   if length(size(H))==3;
      StateVarStructure.H = StateVarStructure.H(:,:,DSGEModel.FirstPeriod+DSGEModel.KalmanFirstObservation-1:DSGEModel.LastPeriod);
   end;
else;
   if length(size(H))==3;
      StateVarStructure.H = StateVarStructure.H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod);
   end;
end;

%
% end of CalculateDSGEStateVariables.m
%
