classdef Monetary < OptimProblem
  properties 
    name          = 'Monetary';
    
    calendar      = [];             % (Optional) calendar associated with the above series
    series        = [];       
    regimes       = [];             % The "estimated" regimes
    
    cond_dist     = 'MarkovChain';  % The (conditional) distribution of the innovations
    method        = 'MLE';          % The estimation method

    wD            = 3;              % Window for U regime
    wU            = 3;              % Window for D regime
    
    %%#####  INTERNAL  ############################################################
    mle  = struct();                % A placeholder for MLE results
  end % end properties
    
  methods 
    %% Construct a Monetary model 
    function self = Monetary(varargin)
     self = self@OptimProblem();
      for no = 1:2:length(varargin)
        setfield(self, varargin{no}, varargin{no+1});
      end      

      % Monetary policy parameters
      self.addParameter('pUU',   0.8632,  [0.00,  1.00]);
      self.addParameter('pSU',   0.0468,  [0.00,  1.00]);   
      self.addParameter('pSD',   0.0647,  [0.00,  1.00]);  
      self.addParameter('pDD',   0.7600,  [0.00,  1.00]);  
      self.addParameter('m0',    3.2,     [   1,  3.99]);  
	  
	  
      self.getRegimes();
    end % end Monetary
    
    %% This function obtains the regimes computed from the reference rate
    function m = getRegimes(self)
      % Lazy copies of wU and wD
      wU = self.wU;
      wD = self.wD;

      T = length(self.series);
      m = NaN(T,1);
      
      % We determine the regime based on Eq. (1) of Begin (2021)
      for dt = 1:T
        Uu = dt < wU || self.series(dt)-self.series(max(1,dt-1)) > 0;
        Uo = (T - dt) < wU || self.series(min(T,dt+1))-self.series(dt) > 0;
        for dU = 2:wU
          Uu = Uu || self.series(dt) - self.series(max(1,dt-dU)) > 0;
          Uo = Uo || self.series(min(T,dt+dU)) - self.series(dt) > 0;
        end
        U = Uu && Uo;

        Du = dt <= wD || self.series(dt)-self.series(dt-1) < 0;
        Do = (T - dt) < wD || self.series(dt+1)-self.series(dt) < 0;
        for dD = 2:wD
          Du = Du || self.series(dt) - self.series(dt-dD) < 0;
          Do = Do || self.series(dt+dU) - self.series(dt) < 0;
        end
        D = Du && Do;

        m(dt) = 2;
        if U
          m(dt) = 1;
        end
        if D
          m(dt) = 3;
        end
        if U && D
          m(dt) = 2;
        end
      end 
      self.regimes = m;
    end % end getRegimes

    %% This functions uses fminsearch to find the MLE parameters
    function [results] = fminsearch(self, lambda, varargin)
      results = fminsearch@OptimProblem(self, lambda, varargin);
      
      switch self.method
        case 'MLE'
          self.mle.params = self.getPV;
          self.mle.out = results;
      end
    end % end fminsearch

    %% This function computes the objective function; the only case considered so far is the MLE
    function S = objective(self, x, varargin)
      switch self.method
        case 'MLE'
          logl = self.getMLEfunction(x);
          S = -sum(logl);
      end
    end % end objective

    %% This function computes the likelihood function for the monetary policy model
    function logl = getMLEfunction(self,p)
      self.setPValues(p);
      pv = self.getPV;
      
      if isempty(self.regimes)
        error('You must filter the regimes first.');
      end
      
      m = [floor(pv.m0);self.regimes];
      
      % We get the right probability from the Markov chain
      T = length(m);
      logl = NaN(T-1,1);
      for dt = 2:T
        if m(dt-1) == 1
          if m(dt) == 1
            logl(dt-1) = log(pv.pUU); 
          elseif m(dt) == 2
            logl(dt-1) = log(1-pv.pUU); 
          else
            logl(dt-1) = log(0); 
          end
        elseif m(dt-1) == 2
          if m(dt) == 1
            logl(dt-1) = log(pv.pSU); 
          elseif m(dt) == 2
            logl(dt-1) = log(1-pv.pSU-pv.pSD); 
          else
            logl(dt-1) = log(pv.pSD); 
          end
        else
          if m(dt) == 1
            logl(dt-1) = log(0);
          elseif m(dt) == 2
            logl(dt-1) = log(1-pv.pDD);
          else
            logl(dt-1) = log(pv.pDD);  
          end
        end
      end
    end % end getMLEfunction
    
    %% This function can generate paths of simulated regimes.
    function [series] = getSimulatedPaths(self,nPaths,nPeriods,initValue)
      pv = self.getPV();
      
      series = NaN(nPaths,nPeriods+1);
      series(:,1) = initValue;
      
      for dt = 2:(nPeriods+1)
        series(:,dt) = (series(:,dt-1) == 1).*randsample(1:3,nPaths,true,[pv.pUU,1-pv.pUU,0])' + ...
                       (series(:,dt-1) == 2).*randsample(1:3,nPaths,true,[pv.pSU,1-pv.pSU-pv.pSD,pv.pSD])' + ...
                       (series(:,dt-1) == 3).*randsample(1:3,nPaths,true,[0,1-pv.pDD,pv.pDD])';
      end
    end % end getSimulatedPaths

    %% This function computes the log prior for the dividend yield model
    function logprior = getlogPrior(self,pv)
      logprior = log((pv.pUU < 1).*(pv.pUU > 0)) + ...
                 log((pv.pSU < 1).*(pv.pSU > 0)) + ...
                 log((pv.pSD < 1).*(pv.pSD > 0)) + ...
                 log((pv.pDD < 1).*(pv.pDD > 0)) + ...
                 log((pv.m0 == 1 | pv.m0 == 2 |  pv.m0 == 3).*1);
    end % end getPrior
	
  end % end methods 

end % end Monetary

