function [Status,Y,YScenario,OriginalShocks,ScenarioShocks,OriginalStates,ScenarioStates,KeepShocks] = DSGEParameterScenariosTheta(DSGEModel,theta,thetaScenario,thetaPositions,ModelParameters,FirstPeriod,LastPeriod,BreakPeriod,CopyFilesToTmpDir,CurrINI)
% DSGEParameterScenariosTheta: Calculate an alternative path for the observed variables under the scenario
%                              parameter values.
%
% USAGE:
%
%       [Status,Y,YScenario,OriginalShocks,ScenarioShocks,OriginalStates,ScenarioStates,KeepShocks] = DSGEParameterScenariosTheta(DSGEModel,theta,thetaScenario,thetaPositions,ModelParameters,FirstPeriod,LastPeriod,BreakPeriod,CopyFilesToTmpDir,CurrINI)
%
% REQUIRED INPUT:  DSGEModel (structure)  whose fields contain all required information
%                            about the model, such as input files, variable and sample
%                            data. Actual observations are not stored in this file.
%
%                  theta (vector) with the original parameter values.
%
%                  thetaScenario (vector) with the scenario values for the 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.
%
%                  FirstPeriod (integer) with the first period in the estimation sample when the training
%                            sample is not taken into account.
%
%                  LastPeriod (integer) with the last period in the estimation sample.
%
%                  BreakPeriod (integer) value between FirstPeriod plus the length of the training sample
%                            and LastPeriod.
%
%                  CopyFilesToTmpDir (boolean) that indicates if certain DSGE model related m-files
%                            should be copied to the tmp directory.
%
%                  CurrINI (structure) with initialization information.
%
% REQUIRED OUTPUT; Status (boolean) which 1 if the calculations were completed without error and
%                            0 otherwise.
%
%                  Y (matrix) with the observed data over the sample for which the scenario is carried out.
%
%                  YScenario (matrix) with the scenario data over the scenario sample.
%
%                  OriginalShocks (matrix) with the smooth estimates of the economic shocks for the original
%                            parameters.
%
%                  ScenarioShocks (matrix) with the smooth estimates of the economic shocks for the scenario
%                            values of the parameters.
%
%                  OriginalStates (matrix) with the smooth estimates of the state variables for the original
%                            parameters.
%
%                  ScenarioStates (matrix) with the state variables under the parameter scenario.
%
%                  KeepShocks (vector) with integers indicating which shocks have a non-zero effect on the
%                            observed variables..
%
%
%                       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: October 26, 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:
%
% * 29-10-2007: Added the "OriginalShocks" and "ScenarioShocks" output arguments.
%
% * 30-10-2007: Updated the documentation.
%
% * 02-11-2007: Added the "KeepShocks" output.
%
% * 05-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.
%
% * 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.
%
% * 07-07-2008: Added code to allow for the case when the measurement matrix H is time-varying.
%
% * 25-07-2008: Added the "BreakPeriod" variable. The scenarios now use the old measurement errors
%               prior the the break period, and the new measurement errors afters the break period.
%               This is consistent with having parameter changes affect the measurement errors.
%
% * 29-07-2008: Added the output variables "OriginalStates" and "ScenarioStates".
%
% * 28-10-2008: Updated the documentation.
%
% * 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
%
Status = 0;
Y = [];
YScenario = [];
OriginalShocks = [];
ScenarioShocks = [];
OriginalStates = [];
ScenarioStates = [];
KeepShocks = [];
%
if CopyFilesToTmpDir==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;
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;
%
% check if the loading of the AiM related data worked without error
%
if isempty(ErrorStr)==0;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% we should now be ready to try to solve the model
%
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 for the original value of the parameters: ' num2str(mcode) ', i.e., "' mcodeStr '" The parameter scenario is therefore aborted.'];
   About(txt,'information','Parameter Scenarios - AiM Solver',200,500,CurrINI);
   return;
end;
if length(size(H))==3;
   if size(H,3)<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(LastPeriod) '). YADA has therefore aborted from the parameter scenario.'];
      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;
%
% now we run the Kalman filter without taking the training sample into account
%
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(:,FirstPeriod:LastPeriod),DSGEModel.X(:,FirstPeriod: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(:,FirstPeriod: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(:,FirstPeriod:LastPeriod),DSGEModel.X(:,FirstPeriod: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(:,FirstPeriod:LastPeriod),[],A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.AllowUnitRoot,1,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(:,FirstPeriod:LastPeriod),DSGEModel.X(:,FirstPeriod:LastPeriod),A,H(:,:,FirstPeriod: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(:,FirstPeriod:LastPeriod),[],A,H(:,:,FirstPeriod: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(:,FirstPeriod:LastPeriod),DSGEModel.X(:,FirstPeriod:LastPeriod),A,H(:,:,FirstPeriod: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(:,FirstPeriod:LastPeriod),[],A,H(:,:,FirstPeriod:LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.AllowUnitRoot,1,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
      end;
   end;
end;
if stat>0;
   txt = ['The Kalman filter provided the return code: ' int2str(stat) '. The state variables are not stationary at the original parameter values. The parameter scenario is therefore aborted.'];
   About(txt,'information','Parameter Scenarios - Kalman Filter',200,500,CurrINI);
   return;
end;
%
% Run the state smoother
%
try;
   if length(size(H))==2;
      [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmoother(DSGEModel.Y(:,FirstPeriod:LastPeriod),Yhat,Ksitt1,Ptt1,H,F,R);
   else;
      [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmootherHt(DSGEModel.Y(:,FirstPeriod:LastPeriod),Yhat,Ksitt1,Ptt1,H(:,:,FirstPeriod:LastPeriod),F,R);
   end;
catch;
   txt = ['The Kalman smoother provided the following error message when using the original parameter values: ' lasterr '. The parameter scenario is therefore aborted.'];
   About(txt,'information','Parameter Scenarios - Kalman Smoother',200,500,CurrINI);
   return;
end;
%
% we will now try to calculate the shocks
%
[S,KeepVar] = RemoveRedundantColumns(B0);
if isempty(S)==1;
   txt = 'The B0 matrix in the state equation on the economic shocks doesn''t have any non-zero columns at the original parameter values. The parameter scenario is therefore aborted.';
   About(txt,'information','Parameter Scenarios - Economic Shocks',200,500,CurrINI);
   return;
end;
KeepShocks = KeepVar;
%
% now we can estimate the economic shocks
%
etatT = inv(S'*S)*S'*[KsitT(:,1) KsitT(:,2:size(KsitT,2))-(F*KsitT(:,1:size(KsitT,2)-1))];
etatT = etatT(:,DSGEModel.KalmanFirstObservation:size(etatT,2));
OriginalStates = KsitT(:,DSGEModel.KalmanFirstObservation:size(KsitT,2));
%
% compute the initial value for the state variables
%
if DSGEModel.KalmanFirstObservation>1;
   KsiT = KsitT(:,DSGEModel.KalmanFirstObservation-1);
else;
   KsiT = zeros(size(KsitT,1),1);
end;
%
% estimate any measurement errors that are available
%
if rank(R)>0;
   if length(size(H))==2;
      if min(size(DSGEModel.X))>=1;
         wtT = DSGEModel.Y(:,FirstPeriod:LastPeriod)-(A'*DSGEModel.X(:,FirstPeriod:LastPeriod))-(H'*KsitT);
      else;
         wtT = DSGEModel.Y(:,FirstPeriod:LastPeriod)-(H'*KsitT);
      end;
   else;
      wtT = zeros(DSGEModel.n,size(etatT,2));
      if min(size(DSGEModel.X))>=1;
         for t=FirstPeriod:LastPeriod;
            wtT(:,t-FirstPeriod+1) = DSGEModel.Y(:,t)-(A'*DSGEModel.X(:,t))-(H(:,:,t)'*KsitT(:,t-FirstPeriod+1));
         end;
      else;
         for t=FirstPeriod:LastPeriod;
            wtT(:,t-FirstPeriod+1) = DSGEModel.Y(:,t)-(H(:,:,t)'*KsitT(:,t-FirstPeriod+1));
         end;
      end;
   end;
   wtT = wtT(:,DSGEModel.KalmanFirstObservation:size(wtT,2));
else;
   wtT = zeros(DSGEModel.n,size(etatT,2));
end;
%
% store the old values in case BreakPeriod is greater than 1
%
Aold = A;
Hold = H;
Rold = R;
Fold = F;
B0old = S;
wtTold = wtT;
etatTold = etatT;
%
% now we try to solve the model using the scenario parameter values
%
[A,H,R,F,B0,mcode] = SolveDSGEModel(DSGEModel,thetaScenario,thetaPositions,ModelParameters,AIMData);
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 for the scenario values of the parameters: ' num2str(mcode) ', i.e., "' mcodeStr '" The parameter scenario is therefore aborted.'];
   About(txt,'information','Parameter Scenarios - AiM Solver',200,500,CurrINI);
   return;
end;
%
% set up the data for the scenario sample
%
Y = DSGEModel.Y(:,FirstPeriod+DSGEModel.KalmanFirstObservation-1:LastPeriod);
if min(size(DSGEModel.X))>=1;
   X = DSGEModel.X(:,FirstPeriod+DSGEModel.KalmanFirstObservation-1:LastPeriod);
end;
%
% remove any zero columns from B0
%
B0 = B0(:,KeepVar);
%
% the economic shocks under the original parameter values
%
OriginalShocks = etatT;
%
% now we will attempt to compute the ScenarioShocks
%
%
% 1. we run the Kalman filter without taking the training sample into account
%
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(:,FirstPeriod:LastPeriod),DSGEModel.X(:,FirstPeriod: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(:,FirstPeriod: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(:,FirstPeriod:LastPeriod),DSGEModel.X(:,FirstPeriod: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(:,FirstPeriod: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(:,FirstPeriod:LastPeriod),DSGEModel.X(:,FirstPeriod:LastPeriod),A,H(:,:,FirstPeriod: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(:,FirstPeriod:LastPeriod),[],A,H(:,:,FirstPeriod: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(:,FirstPeriod:LastPeriod),DSGEModel.X(:,FirstPeriod:LastPeriod),A,H(:,:,FirstPeriod: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(:,FirstPeriod:LastPeriod),[],A,H(:,:,FirstPeriod: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;
if stat>0;
   txt = ['The Kalman filter provided the return code: ' int2str(stat) '. The state variables are not stationary at the scenario parameter values. The parameter scenario is therefore aborted.'];
   About(txt,'information','Parameter Scenarios - Kalman Filter',200,500,CurrINI);
   return;
end;
%
% 2. now we run the Kalman smoother at the scenario values
%
try;
   if length(size(H))==2;
      [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmoother(DSGEModel.Y(:,FirstPeriod:LastPeriod),Yhat,Ksitt1,Ptt1,H,F,R);
   else;
      [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmootherHt(DSGEModel.Y(:,FirstPeriod:LastPeriod),Yhat,Ksitt1,Ptt1,H(:,:,FirstPeriod:LastPeriod),F,R);
   end;
catch;
   txt = ['The Kalman smoother provided the following error message when using the scenario parameter values: ' lasterr '. The parameter scenario is therefore aborted.'];
   About(txt,'information','Parameter Scenarios - Kalman Smoother',200,500,CurrINI);
   return;
end;
[S,KeepVar] = RemoveRedundantColumns(B0);
if size(S)~=size(B0);
   txt = 'The number of non-zero columns of B0 under the parameter scenario values is different from the number of non-zero columns of B0 under the original parameter values. The parameter scenario is therefore aborted.';
   About(txt,'information','Parameter Scenarios - Economic Shocks',200,500,CurrINI);
   return;
end;
%
% now we can estimate the economic shocks
%
etatT = inv(S'*S)*S'*[KsitT(:,1) KsitT(:,2:size(KsitT,2))-(F*KsitT(:,1:size(KsitT,2)-1))];
ScenarioShocks = etatT(:,DSGEModel.KalmanFirstObservation:size(etatT,2));
if BreakPeriod>1;
   ScenarioShocks(:,1:BreakPeriod-1) = etatTold(:,1:BreakPeriod-1);
end;
%
% check if we need to reestimate the measurement errors
%
if rank(R)>0;
   if length(size(H))==2;
      if min(size(DSGEModel.X))>=1;
         wtT = DSGEModel.Y(:,FirstPeriod:LastPeriod)-(A'*DSGEModel.X(:,FirstPeriod:LastPeriod))-(H'*KsitT);
      else;
         wtT = DSGEModel.Y(:,FirstPeriod:LastPeriod)-(H'*KsitT);
      end;
   else;
      wtT = zeros(size(Y));
      if min(size(DSGEModel.X))>=1;
         for t=FirstPeriod:LastPeriod;
            wtT(:,t-FirstPeriod+1) = DSGEModel.Y(:,t)-(A'*DSGEModel.X(:,t))-(H(:,:,t)'*KsitT(:,t-FirstPeriod+1));
         end;
      else;
         for t=FirstPeriod:LastPeriod;
            wtT(:,t-FirstPeriod+1) = DSGEModel.Y(:,t)-(H(:,:,t)'*KsitT(:,t-FirstPeriod+1));
         end;
      end;
   end;
   wtT = wtT(:,DSGEModel.KalmanFirstObservation:size(wtT,2));
else;
   wtT = zeros(size(Y));
end;
%
% now we can calculate the scenario
%
KsiMatrix = zeros(length(KsiT),size(Y,2));
%
if BreakPeriod==1;
   KsiMatrix(:,1) = (F*KsiT)+(B0*etatTold(:,1));
   for i=2:size(etatTold,2);
      KsiMatrix(:,i) = (F*KsiMatrix(:,i-1))+(B0*etatTold(:,i));
   end;
else;
   KsiMatrix(:,1:BreakPeriod-1) = OriginalStates(:,1:BreakPeriod-1);
   for i=BreakPeriod:size(etatTold,2);
      KsiMatrix(:,i) = (F*KsiMatrix(:,i-1))+(B0*etatTold(:,i));
   end;
end;
ScenarioStates = KsiMatrix;
if length(size(H))==2;
   if BreakPeriod==1;
      if min(size(DSGEModel.X))>=1;
         YScenario = (A'*X)+(H'*KsiMatrix)+wtT;
      else;
         YScenario = (H'*KsiMatrix)+wtT;
      end;
   else;
      YScenario = zeros(size(Y));
      if min(size(DSGEModel.X))>=1;
         for t=1:BreakPeriod-1;
            YScenario(:,t) = (Aold'*X(:,t))+(Hold'*KsiMatrix(:,t))+wtTold(:,t);
         end;
         for t=BreakPeriod:size(Y,2);
            YScenario(:,t) = (A'*X(:,t))+(H'*KsiMatrix(:,t))+wtT(:,t);
         end;
      else;
         for t=1:BreakPeriod-1;
            YScenario(:,t) = (Hold'*KsiMatrix(:,t))+wtTold(:,t);
         end;
         for t=BreakPeriod:size(Y,2);
            YScenario(:,t) = (H'*KsiMatrix(:,t))+wtT(:,t);
         end;
      end;
   end;
else;
   YScenario = zeros(size(Y));
   if BreakPeriod==1;
      if min(size(DSGEModel.X))>=1;
         for t=1:size(Y,2);
            YScenario(:,t) = (A'*X(:,t))+(H(:,:,FirstPeriod+DSGEModel.KalmanFirstObservation-2+t)'*KsiMatrix(:,t))+wtT(:,t);
         end;
      else;
         for t=1:size(Y,2);
            YScenario(:,t) = (H(:,:,FirstPeriod+DSGEModel.KalmanFirstObservation-2+t)'*KsiMatrix(:,t))+wtT(:,t);
         end;
      end;
   else;
      if min(size(DSGEModel.X))>=1;
         for t=1:BreakPeriod-1;
            YScenario(:,t) = (Aold'*X(:,t))+(Hold(:,:,FirstPeriod+DSGEModel.KalmanFirstObservation-2+t)'*KsiMatrix(:,t))+wtTold(:,t);
         end;
         for t=BreakPeriod:size(Y,2);
            YScenario(:,t) = (A'*X(:,t))+(H(:,:,FirstPeriod+DSGEModel.KalmanFirstObservation-2+t)'*KsiMatrix(:,t))+wtT(:,t);
         end;
      else;
         for t=1:BreakPeriod-1;
            YScenario(:,t) = (Hold(:,:,FirstPeriod+DSGEModel.KalmanFirstObservation-2+t)'*KsiMatrix(:,t))+wtTold(:,t);
         end;
         for t=BreakPeriod:size(Y,2);
            YScenario(:,t) = (H(:,:,FirstPeriod+DSGEModel.KalmanFirstObservation-2+t)'*KsiMatrix(:,t))+wtT(:,t);
         end;
      end;
   end;
end;
%
Status = 1;

%
% end of DSGEParameterScenariosTheta.m
%
