function [UpdatePriorX,UpdatePriorY,UpdatePriorBounds] = GetUpdatePriorDensity(DSGEModel,CurrINI,controls,CurrentChain,NumChains,NumBurnin,NumIter)
% GetUpdatePriorDensity: Tries to estimate prior densities for the parameters that are a function of
%                        the estimated parameters by drawing from the prior distribution of the latter
%                        parameters.
%
% USAGE:
%
%       [UpdatePriorX,UpdatePriorY,UpdatePriorBounds] = GetUpdatePriorDensity(DSGEModel,CurrINI,controls,CurrentChain,NumChains,NumBurnin,NumIter)
%
% REQUIRED INPUT:  DSGEModel (structure) whose fields contain all required information about the model,
%                            such as input files, variable and sample data.
%
%                  CurrINI (structure) with initialization information.
%
%                  controls (structure), whose fields provide handles to all the UI
%                            controls in YADA's main dialog.
%
%                  CurrentChain (integer) with the number of the currently used MCMC draws.
%
%                  NumChains (integer) with the total number of MCMC chains that we can choose from.
%
%                  NumBurnin (integer) with the number of burn-in draws.
%
%                  NumIter (integer) with the total number of posterior draws that would have
%                            been draws based on the user settings.
%
% REQUIRED OUTPUT: UpdatePriorX (matrix) of dimension NumAdditional x NumValues with the
%                            values on the X-axis for the prior densities of the
%                            NumAdditional additional parameters.
%
%                  UpdatePriorY (matrix) of dimension NumAdditional x NumValues with the
%                            density values for the Y-axis of the prior densities of
%                            the NumAdditional additional parameters.
%
%                  UpdatePriorBounds (matrix) of dimension NumAdditional x 2 with suitable
%                            lower and upper bounds for the additional parameters.
%
%
%                       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 27, 2009.
%                        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:
%
% * 02-03-2009: Renamed the output variables such that the name includes "Prior".
%
% * 09-03-2009: Added quotation marks around paths for the dos copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize
%
UpdatePriorX = [];
UpdatePriorY = [];
UpdatePriorBounds = [];
%
% first we check if the prior densities have already been calculated and saved to file
%
AddFile = [DSGEModel.OutputDirectory '\rwm\RWM-Additional-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(CurrentChain,'%0.0f') '.' num2str(NumChains,'%0.0f') '-' num2str(NumBurnin,'%0.0f') '.' num2str(NumIter,'%0.0f') '.mat'];
if FileExist(AddFile)==1;
   %
   % 
   AddData = load(AddFile);
   if isfield(AddData,'UpdatePriorX')==1;
      UpdatePriorX = AddData.UpdatePriorX;
      UpdatePriorY = AddData.UpdatePriorY;
      UpdatePriorBounds = AddData.UpdatePriorBounds;
      return;
   end;
   %
end;
%
% check if we can load the posterior mode data
%
ModeFile = [DSGEModel.OutputDirectory '\mode\PosteriorMode-' DSGEModel.NameOfModel '.mat'];
if FileExist(ModeFile)==1;
   ModeData = load(ModeFile);
   theta = ModeData.theta;
   thetaDist = ModeData.thetaDist;
   LowerBound = ModeData.LowerBound;
   PriorDist = ModeData.PriorDist;
   ParameterNames = ModeData.ParameterNames;
   thetaPositions = ModeData.thetaPositions;
   ModelParameters = ModeData.ModelParameters;
else;
   [InitStatus,theta,thetaPositions,ModelParameters,thetaDist,PriorDist,ParameterNames,thetaIndex,UniformBounds,LowerBound] = InitializeDSGEModelSimulation(DSGEModel,CurrINI)
   if InitStatus==0;
      return;
   end;
   %
   % check if we have a file with parameters to initialize
   %
   if FileExist(DSGEModel.InitializeParameterFile)==1;
      [stat,msg] = dos(['copy /Y "' DSGEModel.InitializeParameterFile '" "' pwd '\tmp"']);
      if DSGEModel.RunInitializeFirst==1;
         try;
            eval(['ModelParameters = ' GetName(DSGEModel.InitializeParameterFile) '(ModelParameters);']);
         catch;
            return;
         end;
      end;
   end;
end;
%
% copy update parameter file to tmp folder
%
[stat,msg] = dos(['copy /Y "' DSGEModel.UpdateParameterFile '" "' pwd '\tmp"']);
%
% make sure the ModelParameters structure is update with the initial values of the
% estimated parameters
%
ModelParameters = ThetaToModelParameters(ModelParameters,theta,thetaPositions);
%
% run the update parameter function
%
eval(['ModelParameters = ' GetName(DSGEModel.UpdateParameterFile) '(ModelParameters);']);
%
% check if we should run the initialize parameter function
%
if (FileExist(DSGEModel.InitializeParameterFile)==1)&(FileExist(ModeFile)==0)&(DSGEModel.RunInitializeFirst==0);
   try;
      eval(['ModelParameters = ' GetName(DSGEModel.InitializeParameterFile) '(ModelParameters);']);
   catch;
      return;
   end;
end;
%
% the ModelParameter structure now only has values that depend on the initial parameter values.
% determine the values for the additional parameters
%
NumParam = length(theta);
NumAdditional = size(ParameterNames.additional,1);
UpdateInitial = zeros(NumAdditional,1);
for j=1:NumAdditional;
   UpdateInitial(j,1) = getfield(ModelParameters,StringTrim(ParameterNames.additional(j,:)));
end;
%
% try to figure out which estimated parameters affect the update parameters
%
UseParameters = zeros(NumParam,1);
NumBeta = 0;
NumLTN = 0;
NumUniform = 0;
BoundsMatrix = zeros(NumParam,2);
for j=1:NumAdditional;
   for i=1:NumParam;
      if thetaDist(i)==1;
         BoundsMatrix(i,:) = [LowerBound(i,1) Inf];
      elseif thetaDist(i)==2;
         NumBeta = NumBeta+1;
         if size(PriorDist.beta,2)==4;
            BoundsMatrix(i,:) = [PriorDist.beta(NumBeta,3) PriorDist.beta(NumBeta,4)];
         else;
            BoundsMatrix(i,:) = [0 1];
         end;
      elseif thetaDist(i)==3;
         BoundsMatrix(i,:) = [-Inf Inf];
      elseif thetaDist(i)==4;
         BoundsMatrix(i,:) = [LowerBound(i,1) Inf];
      elseif thetaDist(i)==5;
         NumLTN = NumLTN+1;
         BoundsMatrix(i,:) = [PriorDist.beta(NumLTN,3) Inf];
      elseif thetaDist(i)==6;
         NumUniform = NumUniform+1;
         BoundsMatrix(i,:) = [PriorDist.beta(NumUniform,1) PriorDist.beta(NumUniform,2)];
      end;
      thetaNew = theta;
      thetaBounds = sort([0.9*thetaNew(i) 1.1*thetaNew(i)]);
      if thetaBounds(2)<=BoundsMatrix(i,2);
         thetaNew(i) = thetaBounds(2);
      else;
         thetaNew(i) = thetaBounds(1);
      end;
      %
      % update the Model Parameters structure with the thetaNew
      %
      ModelParameters = ThetaToModelParameters(ModelParameters,thetaNew,thetaPositions);
      %
      % run the update parameter function
      %
      eval(['ModelParameters = ' GetName(DSGEModel.UpdateParameterFile) '(ModelParameters);']);
      %
      % check if the initial parameters have changed
      %
      UpdateNew = zeros(NumAdditional,1);
      for j1=1:NumAdditional;
         UpdateNew(j1,1) = getfield(ModelParameters,StringTrim(ParameterNames.additional(j1,:)));
      end;
      if isequal(UpdateNew,UpdateInitial)==0;
         UseParameters(i,1) = i;
      end;
   end;
end;
%
% get only the nonzero values
%
UseParameters = UseParameters(UseParameters>0);
%
% draw from the prior distributions of the parameters to use
%
UseParam = 1;
NumGamma = 0;
NumBeta = 0;
NumNormal = 0;
NumInvGamma = 0;
NumLTN = 0;
NumUniform = 0;
thetaDraws = theta*ones(1,5000);
for i=1:NumParam;
   if thetaDist(i)==1;
      %
      % gamma dist
      %
      NumGamma = NumGamma+1;
      if (UseParameters(UseParam)==i);
         thetaDraws(i,:) = LowerBound(i)+GammaRndFcn(PriorDist.gamma_ab(NumGamma,1),PriorDist.gamma_ab(NumGamma,2),size(thetaDraws,2));
         UseParam = UseParam+1;
      end;
   elseif thetaDist(i)==2;
      %
      % beta dist
      %
      NumBeta = NumBeta+1;
      if size(PriorDist.beta,2)==4;
         %
         % get the lower and upper bounds (if available from PriorDist.)
         %
         c = PriorDist.beta(NumBeta,3);
         d = PriorDist.beta(NumBeta,4);
      else;
         c = 0;
         d = 1;
      end;
      if (UseParameters(UseParam)==i);
         thetaDraws(i,:) = BetaRndFcn(PriorDist.beta_ab(NumBeta,1),PriorDist.beta_ab(NumBeta,2),c,d,size(thetaDraws,2));
         UseParam = UseParam+1;
      end;
   elseif thetaDist(i)==3;
      %
      % normal
      %
      NumNormal = NumNormal+1;
      if (UseParameters(UseParam)==i);
         thetaDraws(i,:) = NormalRndFcn(PriorDist.normal(NumNormal,1),PriorDist.normal(NumNormal,2),size(thetaDraws,2));
         UseParam = UseParam+1;
      end;
   elseif thetaDist(i)==4;
      %
      % inverse gamma
      %
      NumInvGamma = NumInvGamma+1;
      if (UseParameters(UseParam)==i);
         thetaDraws(i,:) = LowerBound(i)+InvGammaRndFcn(PriorDist.invgamma(NumInvGamma,1)-LowerBound(i),PriorDist.invgamma(NumInvGamma,2),size(thetaDraws,2));
         UseParam = UseParam+1;
      end;
   elseif thetaDist(i)==5;
      %
      % left truncated normal
      %
      NumLTN = NumLTN+1;
      if (UseParameters(UseParam)==i);
         thetaDraws(i,:) = LTNormalRndFcn(PriorDist.truncnormal(NumLTN,1),PriorDist.truncnormal(NumLTN,2),PriorDist.truncnormal(NumLTN,3),size(thetaDraws,2));
         UseParam = UseParam+1;
      end;
   elseif thetaDist(i)==6;
      %
      % uniform
      %
      NumUniform = NumUniform+1;
      if (UseParameters(UseParam)==i);
         thetaDraws(i,:) = UniformRndFcn(PriorDist.uniform(NumUniform,1),PriorDist.uniform(NumUniform,2),size(thetaDraws,2));
         UseParam = UseParam+1;
      end;
   end;
   drawnow;
   %
   % check if we're done
   %
   if UseParam>length(UseParameters);
      break;
   end;
end;
%
% we now turn to computing the values of the additional parameters under the draws from the prior
% of the estimated parameters
%
UpdateDraws = zeros(NumAdditional,size(thetaDraws,2));
%
for i=1:size(thetaDraws,2);
   %
   % update ModelParameters
   %
   ModelParameters = ThetaToModelParameters(ModelParameters,thetaDraws(:,i),thetaPositions);
   %
   % run the update parameters function
   %
   eval(['ModelParameters = ' GetName(DSGEModel.UpdateParameterFile) '(ModelParameters);']);
   %
   % collect the update parameter values
   %
   for j=1:NumAdditional;
      UpdateDraws(j,i) = getfield(ModelParameters,StringTrim(ParameterNames.additional(j,:)));
   end;
   drawnow;
end;
%
% determine upper and lower bounds for the update parameters
%
LowerB = floor(min(UpdateDraws'));
UpperB = ceil(max(UpdateDraws'));
if size(LowerB,2)>1
   LowerB = LowerB';
end;
if size(UpperB,2)>1;
   UpperB = UpperB';
end;
UpdatePriorBounds = [LowerB UpperB];
%
% we may now turn to the Kernel density estimation
%
for i=1:NumAdditional;
   [XValues,YValues] = KernelDensityPlot(UpdateDraws(i,:),DSGEModel.KernelDensityEstimator);
   if i==1;
      UpdatePriorX = zeros(NumAdditional,length(XValues));
      UpdatePriorY = UpdatePriorY;
   end;
   UpdatePriorX(i,:) = XValues;
   UpdatePriorY(i,:) = YValues;
   drawnow;
end;
%
% save the results to file so that they can be quickly retrieved
%
save(AddFile,'UpdatePriorX','UpdatePriorY','UpdatePriorBounds','-APPEND');

%
% GetUpdatePriorDensity.m
%
