function [SVEs,status] = DSGECorrelationSimulationTheta(theta,thetaPositions,ModelParameters,NumPaths,EstType,DSGEModel,CurrINI)
% DSGECorrelationSimulationTheta: Estimates observed variable correlations for a DSGE model for a given parameter value.
%                                 The function relies on simulated data and, hence, provides simulation based estimates
%                                 of the sample moments rather than of the population equivalents.
%
% USAGE:
%
%       [SVEs,status] = DSGECorrelationSimulationTheta(theta,thetaPositions,ModelParameters,NumPaths,EstType,DSGEModel,CurrINI)
%
% REQUIRED INPUT:  theta (vector) with the posterior mode estimates of 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.
%
%                  NumPaths (integer) with the number of data paths to simulate per parameter value.
%
%                  EstType (string vector) which determines the type of parameter values that are used. This is
%                                 either "Posterior Mode" or "Initial Values".
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure), with initialization information.
%
% REQUIRED OUTPUT: SVEs (structure) with sample means and quantiles of the mean, covariance, and autocorrelations.
%
% 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
%
%
%                       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: February 29, 2008.
%                        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:
%
% * 03-03-2008: Added the code for computing the mean and the quantiles.
%
% * 01-04-2008: Change the AutoCovHorizon value to "min([SampleFreq*3 T-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.
%
% * 09-07-2008: Added code to allow for a time-varying measurement matrix H.
%
% * 09-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
%
SVEs = [];
%
% 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==2;
      status = 0;
   end;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% now we try to solve the model
%
if nargout==1;
   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;
else;
   ErrorStr = '';
   try;
      [A,H,R,F,B0,status] = 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 status~=1;
      return;
   end;
end;
%
% make sure that only non-zero shocks are kept
%
[B0,KeepShocks] = RemoveRedundantColumns(B0);
n = DSGEModel.n;
k = DSGEModel.k;
[r,q] = size(B0);
%
% 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;
%
% 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);
else;
   Rmeasure = [];
   wSim = 0;
   NumMeasure = 0;
end;
%
% Get the sub-sample
%
[DSGEModel.FirstPeriod,DSGEModel.LastPeriod] = CreateSubSample(DSGEModel);
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 observed variable correlations.'];
      About(ErrorStr,'error','Error - Measurement Equation',200,500,CurrINI);
      return;
   end;
end;
%
% 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);
   XSim = DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod);
else;
   X = [];
   XSim = [];
end;
YTmp = Y;
XTmp = X;
%
% compute the autocorrelations for the actual data
%
T = size(Y,2);
if min(size(X))>0;
   SVEs.Data.Mean = (Y*X')/(X*X');
   %
   % de-mean Y using estimate of the mean from the data rather than A'.
   %
   Y = Y-(SVEs.Data.Mean*X);
else;
   SVEs.Data.Mean = (1/T)*sum(Y,2);
   Y = Y-(SVEs.Data.Mean*ones(1,T));
end;
%
% 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 T-20 18]);
SVEs.AutoCovHorizon = (-AutoCovHorizon:1:AutoCovHorizon);
%
AutoCov = (1/T)*Y*Y';
CorrVec = sqrt(diag(AutoCov));
SVEs.Data.Std = CorrVec;
SVEs.Data.Cov = AutoCov;
SVEs.Data.Corr = zeros((n*(n+1))/2,1+(2*AutoCovHorizon));
SVEs.Data.Corr(:,AutoCovHorizon+1) = vech(AutoCov./(CorrVec*CorrVec'));
%
for i=1:AutoCovHorizon;
   Y = YTmp(:,1+i:size(YTmp,2));
   YL = YTmp(:,1:size(YTmp,2)-i);
   T = size(Y,2);
   if min(size(X))>0;
      X = XTmp(:,1+i:size(YTmp,2));
      XL = XTmp(:,1:size(YTmp,2)-i);
      NewMeanY = (Y*X')/(X*X');
      Y = Y-(NewMeanY*X);
      YL = YL-(NewMeanY*XL);
   else;
      Y = Y-(SVEs.Data.Mean*ones(1,T));
      YL = YL-(SVEs.Data.Mean*ones(1,T));
   end;
   %
   AutoCov = (1/T)*Y*YL';
   CorrVec = sqrt(diag((1/T)*Y*Y'));
   CorrVecL = sqrt(diag((1/T)*YL*YL'));
   Corr = AutoCov./(CorrVec*CorrVecL');
   SVEs.Data.Corr(:,AutoCovHorizon+1+i) = vech(Corr');
   SVEs.Data.Corr(:,AutoCovHorizon+1-i) = vech(Corr);
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;
%
% reset number of observations
%
T = DSGEModel.LastPeriod-DSGEModel.FirstPeriod+1;
Tsum = T-DSGEModel.KalmanFirstObservation+1;
%
% begin the loop
%
abort = '';
MeanEstimationTime = 0;
SVE = [];
%
% check if we should setup a wait handle
%
if DSGEModel.ShowProgress==1;
   ProgressStructure.title = ['Progress for ' num2str(NumPaths,'%0.0f') ' simulation paths'];
   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 = ['Autocorrelation Distribution - ' EstType];
   ProgressStructure.CurrINI = CurrINI;
   WaitHandle = ProgressDLG(0,ProgressStructure);
   set(WaitHandle,'Color',get(0,'defaultuicontrolbackgroundcolor'));
   drawnow;
end;
%
% start the simulation of data
%
if length(size(H))==2;
   for p=1:NumPaths;
      tic;
      %
      % fix the progress dialog is necessary
      %
      if DSGEModel.ShowProgress==1;
         abort = get(WaitHandle,'UserData');
         if strcmp(abort,'cancel')==1;
            StopSimulation = 1;
            break;
         else;
            ProgressDLG([p/NumPaths MeanEstimationTime]);
         end;
      end;
      %
      % draw the economic shocks
      %
      etaSim = randn(q,T);
      %
      % draw the measurement errors
      %
      if isempty(Rmeasure)==0;
         wSim = Rmeasure*randn(NumMeasure,T);
      end;
      %
      % draw an initial state
      %
      StateSim = zeros(r,T+1);
      StateSim(:,1) = Ksi0+(S*sqrt(L)*randn(size(L,2),1));
      %
      % generate the state matrix
      %
      for t=1:T;
         StateSim(:,t+1) = (F*StateSim(:,t))+(B0*etaSim(:,t));
      end;
      drawnow;
      %
      % remove the initial value of the state vector
      %
      StateSim = StateSim(:,2:T+1);
      %
      % create the observed variables
      %
      if k>0;
         YSim = (A'*XSim)+(H'*StateSim)+wSim;
      else;
         YSim = (H'*StateSim)+wSim;
      end;
      %
      % compute the sample moments based on the simulated data
      %
      YTmp = YSim(:,DSGEModel.KalmanFirstObservation:T);
      if k>0;
         XTmp = XSim(:,DSGEModel.KalmanFirstObservation:T);
         MeanY = (YTmp*XTmp')/(XTmp*XTmp');
         YTmp = YTmp-(MeanY*XTmp);
      else;
         MeanY = (1/Tsum)*sum(YTmp,2);
         YTmp = YTmp-MeanY;
      end;
      SVE.CorrData(p).MeanY = MeanY;
      SVE.CorrData(p).StdY = sqrt(diag((1/Tsum)*YTmp*YTmp'));
      %
      % and now for the autocovariances
      %
      SVE.CorrData(p).CorrY = zeros((n*(n+1))/2,1+(2*AutoCovHorizon));
      SVE.CorrData(p).CorrY(:,AutoCovHorizon+1) = vech(((1/Tsum)*YTmp*YTmp')./(SVE.CorrData(p).StdY*SVE.CorrData(p).StdY'));
      for i=1:AutoCovHorizon;
         Yt = YTmp(:,1+i:Tsum);
         YtL = YTmp(:,1:Tsum-i);
         if k>0;
            Xt = XTmp(:,1+i:Tsum);
            Yt = Yt-(((Yt*Xt')/(Xt*Xt'))*Xt);
            XtL = XTmp(:,1:Tsum-i);
            YtL = YtL-(((Yt*Xt')/(Xt*Xt'))*XtL);
         else;
            Yt = Yt-MeanY;
            YtL = YtL-MeanY;
         end;
         Corr = ((1/Tsum)*Yt*YtL')./(SVE.CorrData(p).StdY*SVE.CorrData(p).StdY');
         SVE.CorrData(p).CorrY(:,AutoCovHorizon+1+i) = vech(Corr');
         SVE.CorrData(p).CorrY(:,AutoCovHorizon+1-i) = vech(Corr);         
      end;
      MeanEstimationTime = (((p-1)/p)*MeanEstimationTime)+((1/p)*toc);
      drawnow;
   end;
else;
   for p=1:NumPaths;
      tic;
      %
      % fix the progress dialog is necessary
      %
      if DSGEModel.ShowProgress==1;
         abort = get(WaitHandle,'UserData');
         if strcmp(abort,'cancel')==1;
            StopSimulation = 1;
            break;
         else;
            ProgressDLG([p/NumPaths MeanEstimationTime]);
         end;
      end;
      %
      % draw the economic shocks
      %
      etaSim = randn(q,T);
      %
      % draw the measurement errors
      %
      if isempty(Rmeasure)==0;
         wSim = Rmeasure*randn(NumMeasure,T);
      else;
         wSim = zeros(n,T);
      end;
      %
      % draw an initial state
      %
      StateSim = zeros(r,T+1);
      StateSim(:,1) = Ksi0+(S*sqrt(L)*randn(size(L,2),1));
      %
      % generate the state matrix
      %
      for t=1:T;
         StateSim(:,t+1) = (F*StateSim(:,t))+(B0*etaSim(:,t));
      end;
      drawnow;
      %
      % remove the initial value of the state vector
      %
      StateSim = StateSim(:,2:T+1);
      %
      % create the observed variables
      %
      if k>0;
         for t=DSGEModel.FirstPeriod:DSGEModel.LastPeriod;
            YSim(:,t-DSGEModel.FirstPeriod+1) = (A'*XSim(:,t-DSGEModel.FirstPeriod+1))+(H(:,:,t)'*StateSim(:,t-DSGEModel.FirstPeriod+1))+wSim(:,t-DSGEModel.FirstPeriod+1);
         end;
      else;
         for t=DSGEModel.FirstPeriod:DSGEModel.LastPeriod;
            YSim(:,t-DSGEModel.FirstPeriod+1) = (H(:,:,t)'*StateSim(:,t-DSGEModel.FirstPeriod+1))+wSim(:,t-DSGEModel.FirstPeriod+1);
         end;
      end;
      %
      % compute the sample moments based on the simulated data
      %
      YTmp = YSim(:,DSGEModel.KalmanFirstObservation:T);
      if k>0;
         XTmp = XSim(:,DSGEModel.KalmanFirstObservation:T);
         MeanY = (YTmp*XTmp')/(XTmp*XTmp');
         YTmp = YTmp-(MeanY*XTmp);
      else;
         MeanY = (1/Tsum)*sum(YTmp,2);
         YTmp = YTmp-MeanY;
      end;
      SVE.CorrData(p).MeanY = MeanY;
      SVE.CorrData(p).StdY = sqrt(diag((1/Tsum)*YTmp*YTmp'));
      %
      % and now for the autocovariances
      %
      SVE.CorrData(p).CorrY = zeros((n*(n+1))/2,1+(2*AutoCovHorizon));
      SVE.CorrData(p).CorrY(:,AutoCovHorizon+1) = vech(((1/Tsum)*YTmp*YTmp')./(SVE.CorrData(p).StdY*SVE.CorrData(p).StdY'));
      for i=1:AutoCovHorizon;
         Yt = YTmp(:,1+i:Tsum);
         YtL = YTmp(:,1:Tsum-i);
         if k>0;
            Xt = XTmp(:,1+i:Tsum);
            Yt = Yt-(((Yt*Xt')/(Xt*Xt'))*Xt);
            XtL = XTmp(:,1:Tsum-i);
            YtL = YtL-(((Yt*Xt')/(Xt*Xt'))*XtL);
         else;
            Yt = Yt-MeanY;
            YtL = YtL-MeanY;
         end;
         Corr = ((1/Tsum)*Yt*YtL')./(SVE.CorrData(p).StdY*SVE.CorrData(p).StdY');
         SVE.CorrData(p).CorrY(:,AutoCovHorizon+1+i) = vech(Corr');
         SVE.CorrData(p).CorrY(:,AutoCovHorizon+1-i) = vech(Corr);         
      end;
      MeanEstimationTime = (((p-1)/p)*MeanEstimationTime)+((1/p)*toc);
      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;
end;
%
% now we need to set up the distribution results
%
if strcmp(lower(abort),'cancel')==0;
   NumQuants = length(DSGEModel.Percentiles);
   %
   % initialize model based mean, standard deviation, and autocorrelations
   %
   if k>0;
      SVEs.Mean = zeros(n,k);
   else;
      SVEs.Mean = zeros(n,1);
   end;
   SVEs.Std = zeros(n,1);
   SVEs.Corr = zeros((n*(n+1))/2,1+(2*AutoCovHorizon));
   for i=1:NumQuants;
      if k>0;
         SVEs.Quantiles(i).Mean = zeros(n,k);
      else;
         SVEs.Quantiles(i).Mean = zeros(n,1);
      end;
      SVEs.Quantiles(i).Std = zeros(n,1);
      SVEs.Quantiles(i).Corr = zeros((n*(n+1))/2,1+(2*AutoCovHorizon));
      SVEs.Quantiles(i).percent = DSGEModel.Percentiles(i);
   end;
   %
   % setup a wait dialog
   %
   txt = ['Please wait while YADA computes the distributions for the model based sample means and autocorrelations of the observed variables. Current variable: ' StringTrim(DSGEModel.VariableNames(1,:)) '.'];
   WaitHandle = WaitDLG(txt,'information',['Autocorrelation Distribution - ' EstType],500,200,CurrINI,0);
   WaitControls = get(WaitHandle,'UserData');
   %
   i = 1;
   for i=1:n;
      %
      % setup a wait dialog
      %
      txt = ['Please wait while YADA computes the distributions for the model based sample means and autocorrelations of the observed variables. Current variable: ' StringTrim(DSGEModel.VariableNames(i,:)) '.'];
      set(WaitControls.text,'String',txt);
      drawnow;
      %
      % set up a matrix for the data
      %
      if k>0;
         TempMean = zeros(NumPaths,k);
      else;
         TempMean = zeros(NumPaths,1);
      end;
      TempStd = zeros(NumPaths,1);
      for p=1:NumPaths;
         TempMean(p,:) = SVE.CorrData(p).MeanY(i,:);
         TempStd(p,1) = SVE.CorrData(p).StdY(i,:);
      end;
      drawnow;
      %
      % we can now calculate the mean
      %
      SVEs.Mean(i,:) = (1/NumPaths)*sum(TempMean,1);
      SVEs.Std(i,1) = (1/NumPaths)*sum(TempStd,1);
      %
      % now we can sort these and create the quantiles from the sorted
      %
      for j=1:size(TempMean,2)
         TempMeanj = sort(TempMean(:,j));
         if j==1;
            TempStd = sort(TempStd);
         end;
         for l=1:NumQuants;
            if DSGEModel.Percentiles(l)<50;
               QuantVal = ceil((DSGEModel.Percentiles(l)/100)*NumPaths);
            else;
               QuantVal = floor((DSGEModel.Percentiles(l)/100)*NumPaths);
            end;
            SVEs.Quantiles(l).Mean(i,j) = TempMeanj(QuantVal,1);
            if j==1;
               SVEs.Quantiles(l).Std(i,1) = TempStd(QuantVal,1);
            end;
         end;
         drawnow;
      end;
   end;
   %
   % now we need to deal with the autocorrelations
   %
   SVEs.AutoCovHorizon = (-AutoCovHorizon:1:AutoCovHorizon);
   for j=1:(2*AutoCovHorizon)+1;
      %
      % setup a wait dialog
      %
      txt = ['Please wait while YADA computes the distributions for the model based sample means and autocorrelations of the observed variables. Current autocorrelation horizon: ' int2str(j-(AutoCovHorizon+1)) '.'];
      set(WaitControls.text,'String',txt);
      drawnow;
      %
      TempCorr = zeros(NumPaths,(n*(n+1))/2);
      for p=1:NumPaths;
         TempCorr(p,:) = SVE.CorrData(p).CorrY(:,j)';
      end;
      drawnow;
      %
      % we can now calculate the mean
      %
      SVEs.Corr(:,j) = (1/NumPaths)*sum(TempCorr,1)';
      %
      % sort the data and determine the quantile values
      %
      TempCorr = sort(TempCorr);
      for l=1:NumQuants;
         if DSGEModel.Percentiles(l)<50;
            QuantVal = ceil((DSGEModel.Percentiles(l)/100)*NumPaths);
         else;
            QuantVal = floor((DSGEModel.Percentiles(l)/100)*NumPaths);
         end;
         SVEs.Quantiles(l).Corr(:,j) = TempCorr(QuantVal,:)';
      end;
      drawnow;
   end;
   %
   % close the wait dialog
   %
   set(WaitControls.text,'UserData','done');
   delete(WaitHandle);
   drawnow;
   pause(0.02);
else;
   SVEs = [];
end;

%
% end of DSGECorrelationSimulationTheta.m
%
