classdef ESG
  properties 
    name          = '';

    Monetary      = [];         % A placeholder for the Monetary object
    Inflation     = [];         % A placeholder for the Inflation object    
    ShortRate     = [];         % A placeholder for the ShortRate object    
    Dividends     = [];         % A placeholder for the Dividends object    
    IndexReturns  = [];         % A placeholder for the IndexReturns object    
    ForwardRate   = [];         % A placeholder for the ForwardRate object    
    
    nf     = 100000;            % Size of the "final" Markov chain
    burnin = 100000;            % Burn-in period
    step   = 2;                 % Number of iteration between two samples (thinning)

    beta   = 0.05;              % Adaptation parameter
    d      = 10000;             % Adaptation window
  end % end properties
  
  properties (Access = private)
    n      = NaN;                                                                                                                                                                                            % Actual size of the Markov chain, should be nf x step + burnin + 1
  end
  
  methods
    %% This function creates an ESG object
    function self = ESG(varargin)
      %% Construction of an object
      for no = 1:2:length(varargin)
        self = setfield(self, varargin{no}, varargin{no+1});
      end
      self.n = self.nf * self.step + self.burnin + 1;
    end % end ESG

    %% This function creates a parameter chain using MCMC methods
    function [rechains,acceptance,chains,loglikelihood] = runMarkovChain(self)
      % Initializing chains, acceptance, and loglikelihood variables
      % (this depends on the size of the final parameter set)
      chains.Monetary         = NaN(self.n,5);
      acceptance.Monetary     = NaN(self.n,1);
      loglikelihood.Monetary  = NaN(self.n,1);
      
      % We initilize the chains based on the model selected:
      %   1: Wilkie
      %   2: Extended Wilkie
      %   3: ADG
      %   4: New state-of-the-art ESG
      switch self.Inflation.model_restr
        case 1
          chains.Inflation = NaN(self.n,4);
        case 2
          chains.Inflation = NaN(self.n,4);
        case 3
          chains.Inflation = NaN(self.n,4);
        case 4
          chains.Inflation = NaN(self.n,10);
      end
      acceptance.Inflation    = NaN(self.n,1);
      loglikelihood.Inflation = NaN(self.n,1);
      
      switch self.ShortRate.model_restr
        case 1
          chains.ShortRate = NaN(self.n,5);
        case 2
          chains.ShortRate = NaN(self.n,5);
        case 3
          chains.ShortRate = NaN(self.n,5);
        case 4
          chains.ShortRate = NaN(self.n,11);
      end
      acceptance.ShortRate    = NaN(self.n,1);
      loglikelihood.ShortRate = NaN(self.n,1);
      
      switch self.Dividends.model_restr
        case 1
          chains.Dividends = NaN(self.n,6);
        case 2
          chains.Dividends = NaN(self.n,6);
        case 3
          chains.Dividends = NaN(self.n,6);
        case 4
          chains.Dividends = NaN(self.n,12);
      end
      acceptance.Dividends    = NaN(self.n,1);
      loglikelihood.Dividends = NaN(self.n,1);
      
      switch self.IndexReturns.model_restr
        case 1
          chains.IndexReturns = NaN(self.n,2);
        case 2
          chains.IndexReturns = NaN(self.n,2);
        case 3
          chains.IndexReturns = NaN(self.n,6);
        case 4
          chains.IndexReturns = NaN(self.n,8);
      end
      acceptance.IndexReturns    = NaN(self.n,1);
      loglikelihood.IndexReturns = NaN(self.n,1);
      
      switch self.ForwardRate.model_restr
        case 1
          chains.ForwardRate = NaN(self.n,22);
        case 2
          chains.ForwardRate = NaN(self.n,29);
        case 3
          chains.ForwardRate = NaN(self.n,29);
        case 4
          chains.ForwardRate = NaN(self.n,36);
      end
      acceptance.ForwardRate    = NaN(self.n,1);
      loglikelihood.ForwardRate = NaN(self.n,1);
      
      denMon              = -Inf;
      denInf              = -Inf;
      denShr              = -Inf;
      denDiv              = -Inf;
      denRet              = -Inf;
      denFor              = -Inf;
      
      CovMon              = [];
      CovInf              = [];
      CovShr              = [];
      CovDiv              = [];
      CovRet              = [];
      CovFor              = [];
      
      % Using an MCMC scheme, we generate parameter sets
      textprogressbar(self.n,'Progress: ');
      textprogressbar(1);
      for it = 1:self.n
        % Monetary policy
        [pvalues,accept,denMon]         = self.runMarkovChainMonetary(it,denMon,CovMon);
        chains.Monetary(it,:)           = pvalues;
        acceptance.Monetary(it,1)       = accept;
        loglikelihood.Monetary(it,1)    = denMon;
        
        % Hack: We update the covariance matrix once every 100 iterations
        % to speed up the calculations.
        if mod(it,100) == 0
          CovMon = nancov(chains.Monetary(1:it,1:4));
        end

        % Inflation model
        [pvalues,accept,denInf,zqnew]   = self.runMarkovChainInflation(it,denInf,CovInf);
        chains.Inflation(it,:)          = pvalues;
        acceptance.Inflation(it,1)      = accept;
        loglikelihood.Inflation(it,1)   = denInf;
        if accept || it == 1
          zq = zqnew;
        end
        
        % Hack: We update the covariance matrix once every 100 iterations
        % to speed up the calculations.
        if mod(it,100) == 0
          CovInf = nancov(chains.Inflation(1:it,:));
        end
        
        % Short rate model
        [pvalues,accept,denShr,zrnew]   = self.runMarkovChainShortRate(it,denShr,zq,CovShr);
        chains.ShortRate(it,:)          = pvalues;
        acceptance.ShortRate(it,1)      = accept;
        loglikelihood.ShortRate(it,1)   = denShr;
        if accept || it == 1
          zr = zrnew;
        end
        
        % Hack: We update the covariance matrix once every 100 iterations
        % to speed up the calculations.
        if mod(it,100) == 0
          CovShr = nancov(chains.ShortRate(1:it,:));
        end
        
        % Dividend yield model
        [pvalues,accept,denDiv]         = self.runMarkovChainDividends(it,denDiv,zq,zr,CovDiv);
        chains.Dividends(it,:)          = pvalues;
        acceptance.Dividends(it,1)      = accept;
        loglikelihood.Dividends(it,1)   = denDiv;
        
        % Hack: We update the covariance matrix once every 100 iterations
        % to speed up the calculations.
        if mod(it,100) == 0
          CovDiv = nancov(chains.Dividends(1:it,:));
        end
        
        % Index return model
        [pvalues,accept,denRet]         = self.runMarkovChainIndexReturns(it,denRet,CovRet);
        chains.IndexReturns(it,:)       = pvalues;
        acceptance.IndexReturns(it,1)   = accept;
        loglikelihood.IndexReturns(it,1)= denRet;
        
        % Hack: We update the covariance matrix once every 100 iterations
        % to speed up the calculations.
        if mod(it,100) == 0
          CovRet = nancov(chains.IndexReturns(1:it,:));
        end
        
        % Forward rate model
        [pvalues,accept,denFor]         = self.runMarkovChainForwardRate(it,denFor,CovFor);
        chains.ForwardRate(it,:)        = pvalues;
        acceptance.ForwardRate(it,1)    = accept;
        loglikelihood.ForwardRate(it,1) = denFor;
        
        % Hack: We update the covariance matrix once every 100 iterations
        % to speed up the calculations.
        if mod(it,100) == 0
          CovFor = nancov(chains.ForwardRate(1:it,:));
        end
        
        textprogressbar(it);
      end
      
      % We rearrange the Markov chains for convenience
      [rechains,acceptance,loglikelihood] = self.rearrangeChains(chains,acceptance,loglikelihood);
    end
    
    %% This function generates parameter samples for the monetary policy model
    function [pvnew,accept,num] = runMarkovChainMonetary(self,it,den,covar)
      % FACTOR contains the "well-chosen" factors to be used in the normal
      % mixture
      FACTOR = 1e-2.*ones(1,4);
      
      % We generate uniform rv that will be useful later on
      Uad = unifrnd(0,1); Uac = unifrnd(0,1);

      % Value pvold contains the previous value of the parameter set
      pvold = self.Monetary.getPValues;
      dMom = length(pvold);

      % We generate a new value of the parameter set
      if it < self.d
        pvnew   = mvnrnd(pvold(1:4), FACTOR.^2.*eye(dMom-1)./(dMom-1));
      else
        if Uad > self.beta
          try
            pvnew = mvnrnd(pvold(1:4), (2.38)^2.*covar./(dMom-1));
          catch
            pvnew = mvnrnd(pvold(1:4), FACTOR.^2.*eye(dMom-1)./(dMom-1));
          end
        else
          pvnew = mvnrnd(pvold(1:4), FACTOR.^2.*eye(dMom-1)./(dMom-1));
        end
      end
      
      % We sample an initial regime...
      % Note that Regime 1 is impossible for the initial regime
      pvnew(5)  = randsample(2:3,1); 
      
      % We compute the likelihood and the prior associated with the new 
      % parameter set
      llnum     = sum(self.Monetary.getMLEfunction(pvnew));
      lpriornum = self.Monetary.getlogPrior(self.Monetary.getPV);
      
      % We compute Hastings' ratio and decide whether we keep or not the
      % new parameter value
      if ~isinf(lpriornum) && ~isnan(lpriornum) && ~isinf(llnum) && ~isnan(llnum) && imag(llnum) == 0
        logalpha = min( llnum + lpriornum - den , 0);
      else
        logalpha = -Inf;
      end
      if (log(Uac) < logalpha)
        accept = 1;
        self.Monetary.setPValues(pvnew);
        num = llnum + lpriornum;
      else
        accept = 0;
        self.Monetary.setPValues(pvold);
        num = den;
      end
      pvnew = self.Monetary.getPValues;
    end % end runMarkovChainMonetary

    %% This function generates parameter samples for the inflation model
    function [pvnew,accept,num,zq] = runMarkovChainInflation(self,it,den,covar)
      % FACTOR contains the "well-chosen" factors to be used in the normal
      % mixture
      switch self.Inflation.model_restr
        case 1
          FACTOR = [1e-3,1e-2,1e-6,1e-2];
        case 2
          FACTOR = [1e-3,1e-2,1e-6,1e-2];
        case 3
          FACTOR = [1e-3,1e-2,1e-6,1e-2];
        case 4
          FACTOR = [1e-3,1e-3,1e-3,1e-2,1e-6,1e-2,1e-2,1e1,1e-2,1e-6];
      end

      % We generate uniform rv that will be useful later on
      Uad = unifrnd(0,1); Uac = unifrnd(0,1);
      
      % Value pvold contains the previous value of the parameter set
      pvold = self.Inflation.getPValues;
      dInf = length(pvold);

      % We generate a new value of the parameter set
      if it < self.d
        pvnew     = mvnrnd(pvold, diag(FACTOR.^2)./dInf);
      else
        if Uad > self.beta
          try
            pvnew = mvnrnd(pvold, (2.38)^2.*covar./dInf);
          catch
            pvnew = mvnrnd(pvold, diag(FACTOR.^2)./dInf);
          end
        else
          pvnew   = mvnrnd(pvold, diag(FACTOR.^2)./dInf);
        end
      end
      
      % We compute the likelihood and the prior associated with the new 
      % parameter set
      [llnumc,z,sig2] = self.Inflation.getMLEfunction(pvnew);
      llnum           = sum(llnumc);
      lpriornum       = self.Inflation.getlogPrior(self.Inflation.getPV);
      zq              = z./sqrt(sig2);
      
      % We compute Hastings' ratio and decide whether we keep or not the
      % new parameter value
      if ~isinf(lpriornum) && ~isnan(lpriornum) && ~isinf(llnum) && ~isnan(llnum) && imag(llnum) == 0
        logalpha = min( llnum + lpriornum - den , 0);
      else
        logalpha = -Inf;
      end
      if (log(Uac) < logalpha)
        accept = 1;
        self.Inflation.setPValues(pvnew);
        num = llnum + lpriornum;
      else
        accept = 0;
        self.Inflation.setPValues(pvold);
        num = den;
      end
      pvnew = self.Inflation.getPValues;
    end % end runMarkovChainInflation

    %% This function generates parameter samples for the short rate model
    function [pvnew,accept,num,zr] = runMarkovChainShortRate(self,it,den,zq,covar)
      % FACTOR contains the "well-chosen" factors to be used in the normal
      % mixture
      switch self.ShortRate.model_restr
        case 1
          FACTOR = [1e-2,1e-2,1e-2,1e-6,1e-2];
        case 2
          FACTOR = [1e-2,1e-2,1e-2,1e-6,1e-2];
        case 3
          FACTOR = [1e-2,1e-2,1e-2,1e-6,1e-2];
        case 4
          FACTOR = [1e-2,1e-2,1e-2,1e-2,1e-2,1e-6,1e-2,1e-2,1e1,1e-2,1e-6];
      end
      
      % We generate uniform rv that will be useful later on
      Uad = unifrnd(0,1); Uac = unifrnd(0,1);

      % Value pvold contains the previous value of the parameter set
      pvold = self.ShortRate.getPValues;
      dShr = length(pvold);

      % We generate a new value of the parameter set
      if it < self.d
        pvnew     = mvnrnd(pvold, diag(FACTOR.^2)./dShr);
      else
        if Uad > self.beta
          try
            pvnew = mvnrnd(pvold, (2.38)^2.*covar./dShr);
          catch
            pvnew = mvnrnd(pvold, diag(FACTOR.^2)./dShr);
          end
        else
          pvnew   = mvnrnd(pvold, diag(FACTOR.^2)./dShr);
        end
      end

      % We compute the likelihood and the prior associated with the new 
      % parameter set
      [llnumc,z,sig2] = self.ShortRate.getMLEfunction(pvnew,zq);
      llnum           = sum(llnumc);
      lpriornum       = self.ShortRate.getlogPrior(self.ShortRate.getPV);
      zr              = z./sqrt(sig2);
      
      % We compute Hastings' ratio and decide whether we keep or not the
      % new parameter value
      if ~isinf(lpriornum) && ~isnan(lpriornum) && ~isinf(llnum) && ~isnan(llnum) && imag(llnum) == 0
        logalpha = min( llnum + lpriornum - den , 0);
      else
        logalpha = -Inf;
      end
      if (log(Uac) < logalpha)
        accept = 1;
        self.ShortRate.setPValues(pvnew);
        num = llnum + lpriornum;
      else
        accept = 0;
        self.ShortRate.setPValues(pvold);
        num = den;
      end
      pvnew = self.ShortRate.getPValues;
    end % end runMarkovChainShortRate
    
    %% This function generates parameter samples for the dividend yield model
    function [pvnew,accept,num] = runMarkovChainDividends(self,it,den,zq,zr,covar)
      % FACTOR contains the "well-chosen" factors to be used in the normal
      % mixture
      switch self.Dividends.model_restr
        case 1
          FACTOR = [1e-2,1e-2,1e-2,1e-2,1e-4,1e-2];
        case 2
          FACTOR = [1e-2,1e-2,1e-2,1e-2,1e-4,1e-2];
        case 3
          FACTOR = [1e-2,1e-2,1e-2,1e-2,1e-4,1e-2];
        case 4
          FACTOR = [1e-2,1e-2,1e-2,1e-2,1e-2,1e-2,1e-4,1e-2,1e-2,1e1,1e-2,1e-4];
      end

      % We generate uniform rv that will be useful later on
      Uad = unifrnd(0,1); Uac = unifrnd(0,1);

      % Value pvold contains the previous value of the parameter set
      pvold = self.Dividends.getPValues;
      dDiv = length(pvold);
      
      % We generate a new value of the parameter set
      if it < self.d
        pvnew     = mvnrnd(pvold, diag(FACTOR.^2)./dDiv);
      else
        if Uad > self.beta
          try
            pvnew = mvnrnd(pvold, (2.38)^2.*covar./dDiv);
          catch
            pvnew = mvnrnd(pvold, diag(FACTOR.^2)./dDiv);
          end
        else
          pvnew   = mvnrnd(pvold, diag(FACTOR.^2)./dDiv);
        end
      end

      % We compute the likelihood and the prior associated with the new 
      % parameter set
      llnum           = sum(self.Dividends.getMLEfunction(pvnew,zq,zr,self.ShortRate.params.rho_qr.value));
      lpriornum       = self.Dividends.getlogPrior(self.Dividends.getPV,self.ShortRate.params.rho_qr.value);
      
      % We compute Hastings' ratio and decide whether we keep or not the
      % new parameter value
      if ~isinf(lpriornum) && ~isnan(lpriornum) && ~isinf(llnum) && ~isnan(llnum) && imag(llnum) == 0
        logalpha = min( llnum + lpriornum - den , 0);
      else
        logalpha = -Inf;
      end
      if (log(Uac) < logalpha)
        accept = 1;
        self.Dividends.setPValues(pvnew);
        num = llnum + lpriornum;
      else
        accept = 0;
        self.Dividends.setPValues(pvold);
        num = den;
      end
      pvnew = self.Dividends.getPValues;
    end % end runMarkovChainDividends

    %% This function generates parameter samples for the index return model
    function [pvnew,accept,num] = runMarkovChainIndexReturns(self,it,den,covar)
      % FACTOR contains the "well-chosen" factors to be used in the normal
      % mixture
      switch self.IndexReturns.model_restr
        case 1
          FACTOR = [1e-3,1e-3];
        case 2
          FACTOR = [1e-3,1e-3];
        case 3
          FACTOR = [1e-3.*ones(1,3),1e-3.*ones(1,3)];
        case 4
          FACTOR = [1e-3.*ones(1,3),1e-3,1e-2,1e-2,1e1,1e-2];
      end
      
      % We generate uniform rv that will be useful later on
      Uad = unifrnd(0,1); Uac = unifrnd(0,1);

      % Value pvold contains the previous value of the parameter set
      pvold = self.IndexReturns.getPValues;
      dRet = length(pvold);

      % We generate a new value of the parameter set
      if it < self.d
        pvnew     = mvnrnd(pvold, diag(FACTOR.^2)./dRet);
      else
        if Uad > self.beta
          try
            pvnew = mvnrnd(pvold, (2.38)^2.*covar./dRet);
          catch
            pvnew = mvnrnd(pvold, diag(FACTOR.^2)./dRet);
          end
        else
          pvnew   = mvnrnd(pvold, diag(FACTOR.^2)./dRet);
        end
      end

      % We compute the likelihood and the prior associated with the new 
      % parameter set
      llnum           = sum(self.IndexReturns.getMLEfunction(pvnew));
      lpriornum       = self.IndexReturns.getlogPrior(self.IndexReturns.getPV);
      
      % We compute Hastings' ratio and decide whether we keep or not the
      % new parameter value
      if ~isinf(lpriornum) && ~isnan(lpriornum) && ~isinf(llnum) && ~isnan(llnum) && imag(llnum) == 0
        logalpha = min( llnum + lpriornum - den , 0);
      else
        logalpha = -Inf;
      end
      if (log(Uac) < logalpha)
        accept = 1;
        self.IndexReturns.setPValues(pvnew);
        num = llnum + lpriornum;
      else
        accept = 0;
        self.IndexReturns.setPValues(pvold);
        num = den;
      end
      pvnew = self.IndexReturns.getPValues;
    end % end runMarkovChainIndexReturns
    
    %% This function generates parameter samples for the forward rate model
    function [pvnew,accept,num] = runMarkovChainForwardRate(self,it,den,covar)
      % FACTOR contains the "well-chosen" factors to be used in the normal
      % mixture
      switch self.ForwardRate.model_restr
        case 1
          FACTOR = [1e-3.*ones(1,7),1e-6.*ones(1,7), ...
                    5e-2.*ones(1,2),5e-2.*ones(1,2),5e-6.*ones(1,2),1e-2.*ones(1,2)];
        case 2
          FACTOR = [1e-3.*ones(1,7),1e-6.*ones(1,7),1e-2.*ones(1,7), ...
                    5e-2.*ones(1,2),5e-2.*ones(1,2),5e-6.*ones(1,2),1e-2.*ones(1,2)];
        case 3
          FACTOR = [1e-3.*ones(1,7),1e-6.*ones(1,7),1e-2.*ones(1,7), ...
                    5e-2.*ones(1,2),5e-2.*ones(1,2),5e-6.*ones(1,2),1e-2.*ones(1,2)];
        case 4
          FACTOR = [1e-3.*ones(1,7),1e-6.*ones(1,7),1e-2.*ones(1,7),1e-3.*ones(1,7), ...
                    5e-2.*ones(1,2),5e-2.*ones(1,2),5e-6.*ones(1,2),1e-2.*ones(1,2)];
      end
      
      % We generate uniform rv that will be useful later on
      Uad = unifrnd(0,1); Uac = unifrnd(0,1);

      % Value pvold contains the previous value of the parameter set
      pvold = self.ForwardRate.getPValues;
      dFor = length(pvold);

      % We generate a new value of the parameter set
      if it < self.d
        pvnew     = mvnrnd(pvold, diag(FACTOR.^2)./dFor);
      else
        if Uad > self.beta
          try
            pvnew = mvnrnd(pvold, (2.38)^2.*covar./dFor);
          catch
            pvnew = mvnrnd(pvold, diag(FACTOR.^2)./dFor);
          end
        else
          pvnew   = mvnrnd(pvold, diag(FACTOR.^2)./dFor);
        end
      end

      % We compute the likelihood and the prior associated with the new 
      % parameter set
      llnum           = sum(self.ForwardRate.getMLEfunction(pvnew));
      lpriornum       = self.ForwardRate.getlogPrior(self.ForwardRate.getPV);
      
      % We compute Hastings' ratio and decide whether we keep or not the
      % new parameter value
      if ~isinf(lpriornum) && ~isnan(lpriornum) && ~isinf(llnum) && ~isnan(llnum) && imag(llnum) == 0
        logalpha = min( llnum + lpriornum - den , 0);
      else
        logalpha = -Inf;
      end
      if (log(Uac) < logalpha)
        accept = 1;
        self.ForwardRate.setPValues(pvnew);
        num = llnum + lpriornum;
      else
        accept = 0;
        self.ForwardRate.setPValues(pvold);
        num = den;
      end
      pvnew = self.ForwardRate.getPValues;
    end % end runMarkovChainForwardRate

    %% This function rearranges the Markov chain generated by the function runMarkovChain
    function [rechains,reacceptance,reloglikelihood] = rearrangeChains(self,chains,acceptance,loglikelihood)
      % Monetary policy
      names = fieldnames(self.Monetary.getPV);
      names = names(1:end);
      for dn = 1:length(names)
        rechains.Monetary.(names{dn}) = chains.Monetary((self.burnin+self.step):self.step:end,dn);
      end
      reacceptance.Monetary = acceptance.Monetary((self.burnin+self.step):self.step:end,1);
      reloglikelihood.Monetary = loglikelihood.Monetary((self.burnin+self.step):self.step:end,1);

      % Inflation model
      names = fieldnames(self.Inflation.getPV);
      switch self.Inflation.model_restr
        case 1
          names = names([1,4,5,9]);
        case 2
          names = names([1,4,5,9]);
        case 3
          names = names([1,4,5,9]);
        case 4
          names = names(1:end-1);
      end
      for dn = 1:length(names)
        rechains.Inflation.(names{dn}) = chains.Inflation((self.burnin+self.step):self.step:end,dn);
      end
      reacceptance.Inflation = acceptance.Inflation((self.burnin+self.step):self.step:end,1);
      reloglikelihood.Inflation = loglikelihood.Inflation((self.burnin+self.step):self.step:end,1);

      % Short rate model
      names = fieldnames(self.ShortRate.getPV);
      switch self.ShortRate.model_restr
        case 1
          names = names([1,4,5,6,10]);
        case 2
          names = names([1,4,5,6,10]);
        case 3
          names = names([1,4,5,6,10]);
        case 4
          names = names(1:end-1);
      end
      for dn = 1:length(names)
        rechains.ShortRate.(names{dn}) = chains.ShortRate((self.burnin+self.step):self.step:end,dn);
      end
      reacceptance.ShortRate = acceptance.ShortRate((self.burnin+self.step):self.step:end,1);
      reloglikelihood.ShortRate = loglikelihood.ShortRate((self.burnin+self.step):self.step:end,1);

      % Dividend yield model
      names = fieldnames(self.Dividends.getPV);
      switch self.Dividends.model_restr
        case 1
          names = names([1,4,5,6,7,11]);
        case 2
          names = names([1,4,5,6,7,11]);
        case 3
          names = names([1,4,5,6,7,11]);
        case 4
          names = names(1:end-1);
      end
      for dn = 1:length(names)
        rechains.Dividends.(names{dn}) = chains.Dividends((self.burnin+self.step):self.step:end,dn);
      end
      reacceptance.Dividends = acceptance.Dividends((self.burnin+self.step):self.step:end,1);
      reloglikelihood.Dividends = loglikelihood.Dividends((self.burnin+self.step):self.step:end,1);

      % Index return model
      names = fieldnames(self.IndexReturns.getPV);
      switch self.IndexReturns.model_restr
        case 1
          names = names([1,4]);
        case 2
          names = names([1,4]);
        case 3
          names = names(1:6);
        case 4
          names = names([1,2,3,4,7:10]);
      end
      for dn = 1:length(names)
        rechains.IndexReturns.(names{dn}) = chains.IndexReturns((self.burnin+self.step):self.step:end,dn);
      end
      reacceptance.IndexReturns = acceptance.IndexReturns((self.burnin+self.step):self.step:end,1);
      reloglikelihood.IndexReturns = loglikelihood.IndexReturns((self.burnin+self.step):self.step:end,1);

      % Forward rate model
      names = fieldnames(self.ForwardRate.getPV);
      switch self.ForwardRate.model_restr
        case 1
          names = names([1:14,29:end]);
        case 2
          names = names([1:21,29:end]);
        case 3
          names = names([1:21,29:end]);
        case 4
          names = names(1:end);
      end
      for dn = 1:length(names)
        rechains.ForwardRate.(names{dn}) = chains.ForwardRate((self.burnin+self.step):self.step:end,dn);
      end
      reacceptance.ForwardRate = acceptance.ForwardRate((self.burnin+self.step):self.step:end,1);
      reloglikelihood.ForwardRate = loglikelihood.ForwardRate((self.burnin+self.step):self.step:end,1);
    end % end rearrangeChains
    
  end % end methods
  
end % end ESG