function DoneCalc = DSGEStateVariableEstimation(theta,thetaPostSample,thetaPositions,ModelParameters,SelectedParameters,EstType,TotalDraws,CurrChain,DSGEModel,CurrINI)
% DSGEStateVariableEstimation: Estimates state variables, state shocks, measurement errors and 1-step ahead forecast
%                              for a DSGE  model for draws from the posterior distribution of the structural parameters.
%                              The function also estimates first and second moments of the observed variables and
%                              the state variables using the theoretical formulas. Finally, the function estimates
%                              the first and second moments of the state shocks and the measurement errors treating
%                              these estimates as data. In this case the function provides estimates of the population
%                              moments rather than the sample moments.
%
% USAGE:
%
%       DoneCalc = DSGEStateVariableEstimation(theta,thetaPostSample,thetaPositions,ModelParameters,SelectedParameters,EstType,TotalDraws,CurrChain,DSGEModel,CurrINI)
%
% REQUIRED INPUT:  theta (vector) with the posterior mode estimates of the parameters.
%
%                  thetaPostSample (matrix) of dimension NumDrawsxNumParam with draws from the posterior distribution
%                                  of the original DSGE model parameters.
%
%                  thetaPositions (structure) where thetaPositions(i).parameter is a string vector giving the name
%                                  of the model parameter in position i of theta. The field "distribution" has a string
%                                  vector giving the assumed prior distribution of the parameter.
%
%                  ModelParameters (structure) whose fields have names given by the parameter names in the
%                                 DSGE model.
%
%                  SelectedParameters (vector) of dimension NumParam which indicates which parameters can vary from
%                                 one draw to the next (1) and which parameters are fixed at posterior mode (0).
%
%                  EstType (string vector) that describes what type of estimation is the main estimation.
%
%                  TotalDraws (integer) with the total number of parameter draws that could have been used.
%
%                  CurrChain (integer) with the number for the current Markov chain.
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure), with initialization information.
%
% REQUIRED OUTPUT: DoneCalc (boolean) that is 1 if the calculations were completed and 0 otherwise.
%
%
%                       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: November 21, 2007.
%                        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-11-2007: Added the code for testing if the estimation is run from the "Measurement error" function.
%               In this case, the estimation will now abort if the model doesn't have any measurement
%               errors. The measurement error matrix "R" as well as the "KeepShocks" variables are now also saved
%               to the output file.
%
% * 23-11-2007: Edited a header for an About dialog.
%
% * 27-11-2007: Added the calculations of the covariance and correlation matrices.
%
% * 28-11-2007: Updated the documentation.
%
% * 29-11-2007: Made sure that the "X" matrix is also saved to disk.
%
% * 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.
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 28-02-2008: Updated the documentation.
%
% * 01-04-2008: Changed the AutoCovHorizon value to "min([SampleFreq*3 Tsum-20 18])".
%
% * 16-04-2008: The code now checks if the model has a unique convergent solution also at the posterior mode
%               estimates.
%
% * 23-05-2008: Updated the documentation.
%
% * 29-05-2008: Added a test to check if the model can be solved at the mode.
%
% * 02-07-2008: Added the variable "DSGEModel.AllowUnitRoot" for the Kalman filters.
%
% * 09-07-2008: Added code to allow for the case when the measurement matrix H is time-varying.
%
% * 22-01-2009: Updated the documentation.
%
% * 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
%
DoneCalc = 0;
%
% Prepare variables for selecting the parameters that are updated through the
% posterior draws
%
ScenarioParameters = SelectedParameters .* (1:length(theta));
ScenarioParameters = ScenarioParameters(ScenarioParameters>0);
SelParamStr = strrep(num2str(SelectedParameters),' ','');
thetaMode = theta;
%
% 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;
%
% 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;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% create the output directory
%
stat = MakeDir(DSGEModel.OutputDirectory,'state');
if stat~=1;
   txt = ['YADA was for some reason unable to create the directory "' DSGEModel.OutputDirectory '\state". The computation of ' lower(EstType) ' distributions has therefore been aborted.'];
   About(txt,'information','YADA - Directory Creation Problem',180,500,CurrINI);
   return;
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);
%
% set up Y and X
%
Y = DSGEModel.Y(:,DSGEModel.FirstPeriod+DSGEModel.KalmanFirstObservation-1:DSGEModel.LastPeriod);
if isempty(DSGEModel.X)==0;
   X = DSGEModel.X(:,DSGEModel.FirstPeriod+DSGEModel.KalmanFirstObservation-1:DSGEModel.LastPeriod);
else;
   X = [];
end;
T = DSGEModel.LastPeriod-DSGEModel.FirstPeriod+1;
Tsum = T-DSGEModel.KalmanFirstObservation+1;
%
% initialize index for saving
%
NumDraws = size(thetaPostSample,1);
SaveAfterDraws = min(NumDraws,1000);
NumFiles = ceil(NumDraws/SaveAfterDraws);
%
% solve the model for the posterior mode estimate
%
ErrorStr = '';
try;
   [A,H,R,F,B0,mcode] = SolveDSGEModel(DSGEModel,theta,thetaPositions,ModelParameters,AIMData);
catch;
   ErrorStr = ['YADA caught an error when trying to solve the DSGE model. Message caught is:' lasterr];
end;
if isempty(ErrorStr)==0;
   About(ErrorStr,'information','YADA - Error Message',200,500,CurrINI);
   return;
end;
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;
if length(size(H))==3;
   if size(H,3)<DSGEModel.LastPeriod;
      %
      % 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(DSGEModel.LastPeriod) '). YADA has therefore aborted from the computation of the ' lower(EstType) ' distribution.'];
      About(ErrorStr,'error','Error - Measurement Equation',200,500,CurrINI);
      return;
   end;
end;
%
if (strcmp(lower(EstType),'measurement error')==1)|(strcmp(lower(EstType),'measurement error correlation')==1);
   %
   % check if R has rank greater than 0
   %
   if rank(R)==0;
      txt = ['The DSGE model "' DSGEModel.NameOfModel '" doesn''t have any measurement errors.'];
      About(txt,'information','Measurement Error Distribution',180,500,CurrINI);
      return;
   end;
end;
[B0,KeepShocks] = RemoveRedundantColumns(B0);
%
% set up the autocovariance horizon
%
if strcmp(DSGEModel.DataFrequency,'q')==1;
   SampleFreq = 4;
elseif strcmp(DSGEModel.DataFrequency,'m')==1;
   SampleFreq = 12;
else;
   SampleFreq = 1;
end;
AutoCovHorizon = min([SampleFreq*3 Tsum-20 18]);
%
% set up the date vector
%
DateVector = MakeDateVector(DSGEModel.BeginYear,DSGEModel.BeginPeriod,DSGEModel.EndYear,DSGEModel.EndPeriod,DSGEModel.DataFrequency);
DateVector = DateVector(DSGEModel.FirstPeriod+DSGEModel.KalmanFirstObservation-1:DSGEModel.LastPeriod);
%
% begin the loop
%
abort = '';
LastDraw = 0;
MeanEstimationTime = 0;
ThisSave = 0;
SVE = [];
%
% check if we should setup a wait handle
%
if DSGEModel.ShowProgress==1;
   ProgressStructure.title = ['Progress for ' num2str(NumDraws,'%0.0f') ' parameter draws'];
   ProgressStructure.facecolor = CurrINI.progress_facecolor;
   ProgressStructure.startfacecolor = CurrINI.progress_startfacecolor;
   ProgressStructure.edgecolor = CurrINI.progress_edgecolor;
   ProgressStructure.bgcolor = CurrINI.progress_bgcolor;
   ProgressStructure.stop = 0;
   ProgressStructure.clock = DSGEModel.ShowProgressClock;
   ProgressStructure.label = 'Mean estimation time:';
   %
   ProgressStructure.name = [EstType ' Distribution'];
   ProgressStructure.CurrINI = CurrINI;
   WaitHandle = ProgressDLG(0,ProgressStructure);
   set(WaitHandle,'Color',get(0,'defaultuicontrolbackgroundcolor'));
   drawnow;
else;
   txt = ['Please wait while YADA computes the population based ' lower(EstType) ' distribution. Computations started at: ' StringTrim(datestr(now,14))];
   WaitHandle = WaitDLG(txt,'information',[EstType ' Distribution'],500,150,CurrINI,0);
   WaitControls = get(WaitHandle,'UserData');
end;
%
% start the iterations
%
if length(size(H))==2;
   for it=1:NumDraws;
      LastDraw = LastDraw+1;
      if DSGEModel.ShowProgress==1;
         abort = get(WaitHandle,'UserData');
         if strcmp(abort,'cancel')==1;
            break;
         else;
            ProgressDLG([it/NumDraws MeanEstimationTime]);
         end;
      end;
      %
      % Solve the model for the current parameter vector
      %
      theta(ScenarioParameters) = thetaPostSample(it,ScenarioParameters)';
      tic;
      [A,H,R,F,B0,status] = SolveDSGEModel(DSGEModel,theta,thetaPositions,ModelParameters,AIMData);
      if status~=1;
         %
         % use posterior mode estimate as a backup
         %
         [A,H,R,F,B0,status] = SolveDSGEModel(DSGEModel,thetaMode,thetaPositions,ModelParameters,AIMData);
      end;
      B0 = B0(:,KeepShocks);
      %
      % run the Kalman filter to get the estimate of the state forecasts
      %
      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;
      %
      % calculate the smooth estimates
      %
      [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmoother(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),Yhat,Ksitt1,Ptt1,H,F,R);
      %
      % estimate the economic shocks
      %
      etatT = inv(B0'*B0)*B0'*(KsitT(:,2:T)-(F*KsitT(:,1:T-1)));
      etatt = inv(B0'*B0)*B0'*(Ksitt(:,2:T)-(F*Ksitt(:,1:T-1)));
      %
      % as the first observation we assume that xi0T = 0
      %
      etatT = [inv(B0'*B0)*B0'*KsitT(:,1) etatT];
      etatt = [inv(B0'*B0)*B0'*Ksitt(:,1) etatt];
      if rank(R)>0;
         if min(size(DSGEModel.X))>=1;
            wtT = DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod)-(A'*DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod))-(H'*KsitT);
            wtt = DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod)-(A'*DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod))-(H'*Ksitt);
         else;
            wtT = DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod)-(H'*KsitT);
            wtt = DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod)-(H'*Ksitt);
         end;
         SVE(LastDraw).wtT = wtT(:,DSGEModel.KalmanFirstObservation:T);
         SVE(LastDraw).wtt = wtt(:,DSGEModel.KalmanFirstObservation:T);
         %
         % compute the mean, covariance, and correlation
         %
         SVE(LastDraw).Meanw = (1/Tsum)*sum(SVE(LastDraw).wtT,2);
         wtT = SVE(LastDraw).wtT-(SVE(LastDraw).Meanw*ones(1,Tsum));
         SVE(LastDraw).Covw = (1/Tsum)*wtT*wtT';
         SVE(LastDraw).Corrw = SVE(LastDraw).Covw./(sqrt(diag(SVE(LastDraw).Covw))*sqrt(diag(SVE(LastDraw).Covw))');
      end;
      %
      % make sure we only save the obs we are interested in: starting with DSGEModel.KalmanFirstObservation
      %
      SVE(LastDraw).Yhat = Yhat(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).Ksitt1 = Ksitt1(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).Ksitt = Ksitt(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).KsitT = KsitT(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).etatT = etatT(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).etatt = etatt(:,DSGEModel.KalmanFirstObservation:T);
      %
      % compute the mean, covariance, and correlations for the state variables
      %
      [Sigma,stat] = DoublingAlgorithmLyapunov(F,B0*B0',DSGEModel.DAMaximumIterationsValue*100,10^(-(DSGEModel.DAToleranceValue+1)));
      SVE(LastDraw).StdKsi = sqrt(diag(Sigma));
      SVE(LastDraw).CorrKsi = Sigma./(SVE(LastDraw).StdKsi*SVE(LastDraw).StdKsi');
      %
      % mean, covariance, and correlation for the observed variables
      %
      if min(size(DSGEModel.X))>0;
         SVE(LastDraw).MeanY = A';
      end;
      SVE(LastDraw).StdY = sqrt(diag((H'*Sigma*H)+R));
      SVE(LastDraw).CorrY = zeros(DSGEModel.n,DSGEModel.n,AutoCovHorizon+1);
      SVE(LastDraw).CorrY(:,:,1) = ((H'*Sigma*H)+R)./(SVE(LastDraw).StdY*SVE(LastDraw).StdY');
      Fpower = F;
      for i=1:AutoCovHorizon;
         SVE(LastDraw).CorrY(:,:,1+i) = (H'*Fpower*Sigma*H)./(SVE(LastDraw).StdY*SVE(LastDraw).StdY');
         Fpower = Fpower*F;
      end;
      %
      % compute the mean, covariances, and correlations for the state shocks
      %
      SVE(LastDraw).Meaneta = (1/Tsum)*sum(SVE(LastDraw).etatT,2);
      etatT = SVE(LastDraw).etatT-(SVE(LastDraw).Meaneta*ones(1,Tsum));
      SVE(LastDraw).Coveta = (1/Tsum)*etatT*etatT';
      SVE(LastDraw).Correta = SVE(LastDraw).Coveta./(sqrt(diag(SVE(LastDraw).Coveta))*sqrt(diag(SVE(LastDraw).Coveta))');
      %
      MeanEstimationTime = (((it-1)/it)*MeanEstimationTime)+((1/it)*toc);
      %
      % check if we should save to disk
      %
      if (LastDraw==SaveAfterDraws)|(it==NumDraws);
         ThisSave = ThisSave+1;
         LastDraw = 0;
         %
         % set up the file name and then save
         %
         file = [DSGEModel.OutputDirectory '\state\DSGE-SVE-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' SelParamStr '-' num2str(ThisSave,'%0.0f') '-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '.mat'];
         save(file,'SVE','Y','X','DateVector','R','KeepShocks');
         %
         % restore the data structure
         %
         SVE = [];
      end;
      drawnow;
   end;
else;
   for it=1:NumDraws;
      LastDraw = LastDraw+1;
      if DSGEModel.ShowProgress==1;
         abort = get(WaitHandle,'UserData');
         if strcmp(abort,'cancel')==1;
            break;
         else;
            ProgressDLG([it/NumDraws MeanEstimationTime]);
         end;
      end;
      %
      % Solve the model for the current parameter vector
      %
      theta(ScenarioParameters) = thetaPostSample(it,ScenarioParameters)';
      tic;
      [A,H,R,F,B0,status] = SolveDSGEModel(DSGEModel,theta,thetaPositions,ModelParameters,AIMData);
      if status~=1;
         %
         % use posterior mode estimate as a backup
         %
         [A,H,R,F,B0,status] = SolveDSGEModel(DSGEModel,thetaMode,thetaPositions,ModelParameters,AIMData);
      end;
      B0 = B0(:,KeepShocks);
      %
      % run the Kalman filter to get the estimate of the state forecasts
      %
      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;
      %
      % calculate the smooth estimates
      %
      [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmootherHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),Yhat,Ksitt1,Ptt1,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,R);
      %
      % estimate the economic shocks
      %
      etatT = inv(B0'*B0)*B0'*(KsitT(:,2:T)-(F*KsitT(:,1:T-1)));
      etatt = inv(B0'*B0)*B0'*(Ksitt(:,2:T)-(F*Ksitt(:,1:T-1)));
      %
      % as the first observation we assume that xi0T = 0
      %
      etatT = [inv(B0'*B0)*B0'*KsitT(:,1) etatT];
      etatt = [inv(B0'*B0)*B0'*Ksitt(:,1) etatt];
      if rank(R)>0;
         wtT = zeros(DSGEModel.n,size(KsitT,2));
         if min(size(DSGEModel.X))>=1;
            for t=DSGEModel.FirstPeriod:DSGEModel.LastPeriod;
               wtT(:,t-DSGEModel.FirstPeriod+1) = DSGEModel.Y(:,t)-(A'*DSGEModel.X(:,t))-(H(:,:,t)'*KsitT(:,t-DSGEModel.FirstPeriod+1));
               wtt(:,t-DSGEModel.FirstPeriod+1) = DSGEModel.Y(:,t)-(A'*DSGEModel.X(:,t))-(H(:,:,t)'*Ksitt(:,t-DSGEModel.FirstPeriod+1));
            end;
         else;
            for t=DSGEModel.FirstPeriod:DSGEModel.LastPeriod;
               wtT(:,t-DSGEModel.FirstPeriod+1) = DSGEModel.Y(:,t)-(H(:,:,t)'*KsitT(:,t-DSGEModel.FirstPeriod+1));
               wtt(:,t-DSGEModel.FirstPeriod+1) = DSGEModel.Y(:,t)-(H(:,:,t)'*Ksitt(:,t-DSGEModel.FirstPeriod+1));
            end;
         end;
         SVE(LastDraw).wtT = wtT(:,DSGEModel.KalmanFirstObservation:T);
         SVE(LastDraw).wtt = wtt(:,DSGEModel.KalmanFirstObservation:T);
         %
         % compute the mean, covariance, and correlation
         %
         SVE(LastDraw).Meanw = (1/Tsum)*sum(SVE(LastDraw).wtT,2);
         wtT = SVE(LastDraw).wtT-(SVE(LastDraw).Meanw*ones(1,Tsum));
         SVE(LastDraw).Covw = (1/Tsum)*wtT*wtT';
         SVE(LastDraw).Corrw = SVE(LastDraw).Covw./(sqrt(diag(SVE(LastDraw).Covw))*sqrt(diag(SVE(LastDraw).Covw))');
      end;
      %
      % make sure we only save the obs we are interested in: starting with DSGEModel.KalmanFirstObservation
      %
      SVE(LastDraw).Yhat = Yhat(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).Ksitt1 = Ksitt1(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).Ksitt = Ksitt(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).KsitT = KsitT(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).etatT = etatT(:,DSGEModel.KalmanFirstObservation:T);
      SVE(LastDraw).etatt = etatt(:,DSGEModel.KalmanFirstObservation:T);
      %
      % compute the mean, covariance, and correlations for the state variables
      %
      [Sigma,stat] = DoublingAlgorithmLyapunov(F,B0*B0',DSGEModel.DAMaximumIterationsValue*100,10^(-(DSGEModel.DAToleranceValue+1)));
      SVE(LastDraw).StdKsi = sqrt(diag(Sigma));
      SVE(LastDraw).CorrKsi = Sigma./(SVE(LastDraw).StdKsi*SVE(LastDraw).StdKsi');
      %
      % mean, covariance, and correlation for the observed variables
      %
      if min(size(DSGEModel.X))>0;
         SVE(LastDraw).MeanY = A';
      end;
      CovY = zeros(DSGEModel.n,DSGEModel.n);
      for t=DSGEModel.FirstPeriod+DSGEModel.KalmanFirstObservation-1:DSGEModel.LastPeriod;
         CovY = CovY+((1/(DSGEModel.LastPeriod+1-(DSGEModel.FirstPeriod+DSGEModel.KalmanFirstObservation-1)))*H(:,:,t)'*Sigma*H(:,:,t));
      end;
      SVE(LastDraw).StdY = sqrt(diag(CovY+R));
      SVE(LastDraw).CorrY = zeros(DSGEModel.n,DSGEModel.n,AutoCovHorizon+1);
      SVE(LastDraw).CorrY(:,:,1) = (CovY+R)./(SVE(LastDraw).StdY*SVE(LastDraw).StdY');
      Fpower = F;
      for i=1:AutoCovHorizon;
         ACovY = zeros(DSGEModel.n,DSGEModel.n);
         for t=DSGEModel.FirstPeriod+DSGEModel.KalmanFirstObservation-1:DSGEModel.LastPeriod-i;
            ACovY = ACovY+((1/(DSGEModel.LastPeriod+1-i-(DSGEModel.FirstPeriod+DSGEModel.KalmanFirstObservation-1)))*(H(:,:,t+i)'*Fpower*Sigma*H(:,:,t)));
         end;
         SVE(LastDraw).CorrY(:,:,1+i) = ACovY./(SVE(LastDraw).StdY*SVE(LastDraw).StdY');
         Fpower = Fpower*F;
      end;
      %
      % compute the mean, covariances, and correlations for the state shocks
      %
      SVE(LastDraw).Meaneta = (1/Tsum)*sum(SVE(LastDraw).etatT,2);
      etatT = SVE(LastDraw).etatT-(SVE(LastDraw).Meaneta*ones(1,Tsum));
      SVE(LastDraw).Coveta = (1/Tsum)*etatT*etatT';
      SVE(LastDraw).Correta = SVE(LastDraw).Coveta./(sqrt(diag(SVE(LastDraw).Coveta))*sqrt(diag(SVE(LastDraw).Coveta))');
      %
      MeanEstimationTime = (((it-1)/it)*MeanEstimationTime)+((1/it)*toc);
      %
      % check if we should save to disk
      %
      if (LastDraw==SaveAfterDraws)|(it==NumDraws);
         ThisSave = ThisSave+1;
         LastDraw = 0;
         %
         % set up the file name and then save
         %
         file = [DSGEModel.OutputDirectory '\state\DSGE-SVE-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' SelParamStr '-' num2str(ThisSave,'%0.0f') '-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '.mat'];
         save(file,'SVE','Y','X','DateVector','R','KeepShocks');
         %
         % restore the data structure
         %
         SVE = [];
      end;
      drawnow;
   end;
end;
%
% close the wait dialog
%
if DSGEModel.ShowProgress==1;
   if ishandle(WaitHandle)==1;
      set(WaitHandle,'UserData','done');
      close(WaitHandle);
      drawnow;
      pause(0.02);
   end;
else;
   set(WaitControls.text,'UserData','done');
   delete(WaitHandle);
   drawnow;
   pause(0.02);
end;
%
% check if we have computed all output
%
if strcmp(abort,'cancel')==0;
   DoneCalc = 1;
end;

%
% end of DSGEStateVariableEstimation.m
%
