function [SimulatedData,KeepShocks] = DSGESimulateVariables(A,H,R,F,B0,T,DSGEModel,CurrINI)
% DSGESimulateVariables: Simulates a sample of data for the observed variables, the state variables, shocks, and errors
%
% USAGE:
%
%       [SimulatedData,KeepShocks] = DSGESimulateVariables(A,H,R,F,B0,NumObs(SelectedOption),DSGEModel,CurrINI)
%
% REQUIRED INPUT:  A (matrix) of dimension kxn with the parameters on the exogenous variables in
%                      the measurement equation.
%
%                  H (matrix) of dimension rxn with parameters on the state variabels in the
%                      measurement equation.
%
%                  R (matrix) of dimension nxn with the covariances of the measurement errors
%                      in the measurement equation.
%
%                  F (matrix) of dimension rxr with state transitions parameters for the state
%                      equation.
%
%                  B0 (matrix) of dimension rxq with parameters on the economic shocks in the
%                      state equation.
%
%                  T (integer) with the number of observations to simulate.
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure), with initialization information.
%
% REQUIRED OUTPUT: SimulatedData (structure) with fields Y (observed variables), ksi (state variables),
%                      eta (economic shocks), w (measurement errors), and X (exogenous variables). Also,
%                      the output structure includes the field Parameters, with 7 subfields given by,
%                      A, H, R, F, B0, ksi0, and P0. The last two contain the population mean and
%                      covariance for the state variables.
%
%                  KeepShocks (vector) with indices for the economic shocks that are used for the simulated data.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: July 11, 2008.
%                        This version: December 18, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:
%
% * 15-7-2008:  Fixed a bug when setting up a 3D H-matrix.
%
% * 18-12-2008: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize
%
SimulatedData = [];
%
% create deterministic variabels
%
if isempty(DSGEModel.X)==0;
   if T<=DSGEModel.T;
      X = DSGEModel.X;
   else;
      %
      % use the last observations
      %
      if strcmp(lower(DSGEModel.DataFrequency),'q')==1;
         %
         % get the last four observations
         %
         LastX = DSGEModel.X(:,DSGEModel.T-3:DSGEModel.T);
         NumYears = ceil((T-DSGEModel.T)/4);
         NewX = kron(ones(1,NumYears),LastX);
         X = [DSGEModel.X NewX(:,1:T-DSGEModel.T)];
      elseif strcmp(lower(DSGEModel.DataFrequency),'m')==1;
         %
         % get the last twelve observations
         %
         LastX = DSGEModel.X(:,DSGEModel.T-11:DSGEModel.T);
         NumYears = ceil((T-DSGEModel.T)/12);
         NewX = kron(ones(1,NumYears),LastX);
         X = [DSGEModel.X NewX(:,1:T-DSGEModel.T)];
      else;
         X = [DSGEModel.X DSGEModel.X(:,DSGEModel.T)*ones(1,T-DSGEModel.T)];
      end;
   end;
else;
   X = [];
end;
%
% make sure that only the non-zero economic shocks are used
%
[B0,KeepShocks] = RemoveRedundantColumns(B0);
n = DSGEModel.n;
k = DSGEModel.k;
[r,q] = size(B0);
%
% fix the measurement matrix H is time-varying
%
if length(size(H))==3;
   if T<=DSGEModel.T;
      H = H(:,:,1:T);
   else;
      Hold = H;
      H = zeros(r,n,T);
      H(:,:,1:DSGEModel.T) = Hold(:,:,1:DSGEModel.T);
      %
      % use the last observed H value for all future values
      %
      for t=DSGEModel.T+1:T;
         if size(Hold,3)>=t;
            H(:,:,t) = Hold(:,:,t);
         else;
            H(:,:,t) = Hold(:,:,size(Hold,3));
         end;
      end;
   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;
%
% initialize the random number generator
%
if DSGEModel.RandomNumberValue==1;
   randn('state',0);
   rand('state',0);
else;
   randn('state',sum(100*clock));
   rand('state',sum(100*clock));
end;
%
% initial values for the mean and covariance matrix of the distribution of the state variables
%
if isempty(DSGEModel.UnitRootStates)==1;
   [P0,stat] = DoublingAlgorithmLyapunov(F,B0*B0',DSGEModel.DAMaximumIterationsValue*100,10^(-(DSGEModel.DAToleranceValue+1)));
else;
   StationaryPos = GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates);
   P0 = DSGEModel.StateCovConstValue*eye(r,r);
   F1 = F(StationaryPos,StationaryPos);
   Q1 = B0(StationaryPos,:)*B0(StationaryPos,:)';
   [Pstat,stat] = DoublingAlgorithmLyapunov(F1,Q1,DSGEModel.DAMaximumIterationsValue*100,10^(-(DSGEModel.DAToleranceValue+1)));
   P0(StationaryPos,StationaryPos) = Pstat;
end;
ksi0 = DSGEModel.InitialStateVector;
[S,L] = SchurDecomposition(P0);
%
% deal with measurement errors
%
if rank(R)>0;
   [Rm,Lm] = SchurDecomposition(R);
   Rmeasure = Rm*sqrt(Lm);
   NumMeasure = size(Lm,1);
   w = Rmeasure*randn(NumMeasure,T);
else;
   Rmeasure = [];
   w = zeros(n,T);
   NumMeasure = 0;
end;
%
% draw the economic shocks
%
eta = randn(q,T);
shocks = B0*eta;
%
% simulate the state variables
%
ksi = zeros(r,T+1);
ksi(:,1) =  ksi0+(S*sqrt(L)*randn(size(L,2),1));
for t=1:T;
   ksi(:,t+1) = F*ksi(:,t)+shocks(:,t);
end;
ksi = ksi(:,2:T+1);
if length(size(H))==2;
   if isempty(X)==0;
      Y = (A'*X)+(H'*ksi)+w;
   else;
      Y = (H'*ksi)+w;
   end;
else;
   Y = zeros(n,T);
   if isempty(X)==0;
      for t=1:T;
         Y(:,t) = (A'*X(:,t))+(H(:,:,t)'*ksi(:,t))+w(:,t);
      end;
   else;
      for t=1:T;
         Y(:,t) = (H(:,:,t)'*ksi(:,t))+w(:,t);
      end;
   end;
end;
%
% setup the output
%
SimulatedData.Y = Y;
SimulatedData.X = X;
SimulatedData.ksi = ksi;
SimulatedData.w = w;
SimulatedData.eta = eta;
SimulatedData.Parameters.A = A;
SimulatedData.Parameters.H = H;
SimulatedData.Parameters.R = R;
SimulatedData.Parameters.F = F;
SimulatedData.Parameters.B0 = B0;
SimulatedData.Parameters.ksi0 = ksi0;
SimulatedData.Parameters.P0 = P0;

%
% end of DSGESimulateVariables.m
%
