classdef Dividends < OptimProblem
  properties 
    name          = 'Dividends';
    
    calendar      = [];         % (Optional) calendar associated with the above series
    series        = [];       

    Monetary      = [];         % A placeholder for the Monetary object
    Inflation     = [];         % A placeholder for the Inflation object    
    ShortRate     = [];         % A placeholder for the ShortRate object    

    cond_dist     = 'Normal';   % The (conditional) distribution of the innovations. 
    method        = 'MLE';      % The estimation method
    
    model_restr   =  1;         % 1 for Wilkie-like, 
                                % 2 for Extended Wilkie-like, 
                                % 3 for ADG-like, 
                                % 4 for full new ESG
    
    %%#####  INTERNAL  ############################################################
    mle  = struct();            % A placeholder for MLE results
  end % end properties
    
  methods 
    %% Construct a Inflation model 
    function self = Dividends(varargin)
     self = self@OptimProblem();
      for no = 1:2:length(varargin)
        setfield(self, varargin{no}, varargin{no+1});
      end      

      % Inflation parameters
      self.addParameter('mu_d1',           -3.6264,       [-5.000, 0.000]);
      self.addParameter('mu_d2',           -4.0409,       [-5.000, 0.000]);
      self.addParameter('mu_d3',           -3.7146,       [-5.000, 0.000]);
      self.addParameter('a_d',              0.9850,       [-1.000, 1.000]);
      self.addParameter('rho_qd',          -0.1042,       [-1.000, 1.000]);
      self.addParameter('rho_rd',          -0.0537,       [-1.000, 1.000]);
      self.addParameter('sig2_d',           0.0029,       [ 1e-9,  0.005]);
      self.addParameter('alpha_d',          0.1266,       [ 0.00, 1.0000]);
      self.addParameter('alphabetagamma_d', 0.9796,       [ 0.00, 0.9950]);
      self.addParameter('gamma_d',          0.0000,       [  -10,     10]);
      self.addParameter('d0',               0.0578,       [ 0.00, 0.1000]);
      self.addParameter('sig2_d_init',      0.0022,       [ 0.00,   0.05]);

      % Depending on the model selected, we turn on or off some parameters
      if self.model_restr == 1
        self.params.mu_d2.fixed             = true;
        self.params.mu_d3.fixed             = true;
        self.params.alpha_d.fixed           = true;
        self.params.alphabetagamma_d.fixed  = true;
        self.params.gamma_d.fixed           = true;
        self.params.sig2_d_init.fixed       = true;

        self.params.alpha_d.value           = 0;
        self.params.alphabetagamma_d.value  = 0;
        self.params.gamma_d.value           = 0;
        self.params.sig2_d_init.value       = self.params.sig2_d.value;
        self.params.mu_d2.value             = self.params.mu_d1.value;
        self.params.mu_d3.value             = self.params.mu_d1.value;
      elseif self.model_restr == 2
        self.params.mu_d2.fixed             = true;
        self.params.mu_d3.fixed             = true;
        self.params.alpha_d.fixed           = true;
        self.params.alphabetagamma_d.fixed  = true;
        self.params.gamma_d.fixed           = true;
        self.params.sig2_d_init.fixed       = true;

        self.params.alpha_d.value           = 0;
        self.params.alphabetagamma_d.value  = 0;
        self.params.gamma_d.value           = 0;
        self.params.sig2_d_init.value       = self.params.sig2_d.value;
        self.params.mu_d2.value             = self.params.mu_d1.value;
        self.params.mu_d3.value             = self.params.mu_d1.value;
      elseif self.model_restr == 3
        self.params.mu_d2.fixed             = true;
        self.params.mu_d3.fixed             = true;
        self.params.alpha_d.fixed           = true;
        self.params.alphabetagamma_d.fixed  = true;
        self.params.gamma_d.fixed           = true;
        self.params.sig2_d_init.fixed       = true;

        self.params.alpha_d.value           = 0;
        self.params.alphabetagamma_d.value  = 0;
        self.params.gamma_d.value           = 0;
        self.params.sig2_d_init.value       = self.params.sig2_d.value;
        self.params.mu_d2.value             = self.params.mu_d1.value;
        self.params.mu_d3.value             = self.params.mu_d1.value;
      end
    end % end Dividends
    
    %% This function returns the parameters
    function pv = getPV(self)
      pv        = getPV@OptimProblem(self);
      
      % Again, depending on the model selected, some parameters are turned
      % on or off.
      if self.model_restr == 1
        self.params.sig2_d_init.value = self.params.sig2_d.value;
        self.params.mu_d2.value       = self.params.mu_d1.value;
        self.params.mu_d3.value       = self.params.mu_d1.value;
        
        pv.sig2_d_init                = self.params.sig2_d.value;
        pv.mu_d2                      = self.params.mu_d1.value;
        pv.mu_d3                      = self.params.mu_d1.value;
      elseif self.model_restr == 2
        self.params.sig2_d_init.value = self.params.sig2_d.value;
        self.params.mu_d2.value       = self.params.mu_d1.value;
        self.params.mu_d3.value       = self.params.mu_d1.value;
        
        pv.sig2_d_init                = self.params.sig2_d.value;
        pv.mu_d2                      = self.params.mu_d1.value;
        pv.mu_d3                      = self.params.mu_d1.value;
      elseif self.model_restr == 3
        self.params.sig2_d_init.value = self.params.sig2_d.value;
        self.params.mu_d2.value       = self.params.mu_d1.value;
        self.params.mu_d3.value       = self.params.mu_d1.value;
        
        pv.sig2_d_init                = self.params.sig2_d.value;
        pv.mu_d2                      = self.params.mu_d1.value;
        pv.mu_d3                      = self.params.mu_d1.value;
      end
      
      pv.beta_d = pv.alphabetagamma_d - pv.alpha_d*(1+pv.gamma_d^2);
    end % end getPV

    %% This functions uses fminsearch to find the MLE parameters
    function [results] = fminsearch(self, lambda, varargin)
      results = fminsearch@OptimProblem(self, lambda, varargin{:});     
      [~,z,sig2] = self.objective(self.getPValues, varargin{:});
      
      switch self.method
        case 'MLE'
          self.mle.params = self.getPV;
          self.mle.z      = z;
          self.mle.sig2   = sig2;
          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,z,sig2] = objective(self, x, varargin)
      switch self.method
        case 'MLE'
          [logl,z,sig2] = self.getMLEfunction(x,varargin{:});
          S = -sum(logl);
      end
    end % end objective

    %% This function computes the likelihood function for the dividend yield model
    function [logl,z,sig2] = getMLEfunction(self,p,zq,zr,rho_qr)
      % We extract the current value of the parameters
      self.setPValues(p);
      pv = self.getPV;
      pv.rho_qr = rho_qr;

      % We try to compute the Cholesky decomposition of our parameters; if
      % it does not exist, we return -Inf as the log-likelihood to discard
      % this set of parameters.
      try
        C = chol([1,pv.rho_qr,pv.rho_qd;pv.rho_qr,1,pv.rho_rd;pv.rho_qd,pv.rho_rd,1]);
      catch
        logl = -Inf;
        z    = [];
        sig2 = [];
        return;
      end
      
      % We extract the number of observations and the monetary policy chain
      T = length(self.series) + 1;
      m = self.Monetary.regimes; m = [floor(self.Monetary.params.m0.value);m];
      
      z    = zeros(T,1);
      y    = log([pv.d0;self.series]);
      sig2 = zeros(T+1,1);

      mus = [pv.mu_d1,pv.mu_d2,pv.mu_d3];

      % For each observation, we compute the log-likelihood contribution at
      % time t
      sig2(2) = pv.sig2_d_init;
      for dt = 2:T
        z(dt) = y(dt) - mus(m(dt)) - pv.a_d*(y(dt-1) - mus(m(dt)));
        sig2(dt+1) = pv.sig2_d + pv.beta_d*(sig2(dt)-pv.sig2_d) + pv.alpha_d*( (z(dt) - pv.gamma_d.*sqrt(sig2(dt))).^2 - (1+pv.gamma_d.^2)*pv.sig2_d );
      end
      sig2 = sig2(2:end-1);
      z    = z(2:end);
      
      z1   = zq(:);
      z2   = (zr(:) - C(1,2).*zq(:))./C(2,2);
      
      logl = log( normpdf( z - C(1,3).*sqrt(sig2).*z1 - C(2,3).*sqrt(sig2).*z2, 0, sqrt(sig2)*C(3,3) ));
      if imag(sum(logl)) ~= 0
        logl = NaN(size(logl));
      end
    end % end getMLEfunction
    
    %% This function computes the log prior for the dividend yield model
    function logprior = getlogPrior(self,pv,rho_qr)
      pv.rho_qr = rho_qr;
      
      % We need to compute the determinant of our correlation coefficients
      % to obtain the LKJ prior distribution (proportional to the
      % determinant of the correlation matrix)
      S        = [1,pv.rho_qr,pv.rho_qd; ...
                  pv.rho_qr,1,pv.rho_rd; ...
                  pv.rho_qd,pv.rho_rd,1];
      detS     = det(S);
      
      if detS < 0
        logprior = -Inf;
      else
        logprior =  log(normpdf(exp(pv.mu_d1),0,1)) + ...
                    log(normpdf(exp(pv.mu_d2),0,1)) + ...
                    log(normpdf(exp(pv.mu_d3),0,1)) + ...
                    log((pv.a_d <= 0.995).*(pv.a_d >= -0.995)) + ...                  % Hack: Otherwise, we might have some stationarity issues...
                    log(detS) + ...                                                   % Proportional to the LKJ distribution
                    log(exppdf(pv.sig2_d,1)) + ...
                    log((pv.alpha_d <= 0.995).*(pv.alpha_d >= 0)) + ...               % Hack: Otherwise, we might have some stationarity issues...
                    log((pv.alphabetagamma_d <= 0.995).*(pv.alphabetagamma_d >= 0)) + ...       % Hack: Otherwise, we might have some stationarity issues...
                    log((pv.beta_d <= 0.995).*(pv.beta_d >= 0)) + ...    			  % Hack: Otherwise, we might have some stationarity issues...
                    log(normpdf(pv.gamma_d,0,10)) + ...
                    log(normpdf(pv.d0,0,1)) + ...
                    log(exppdf(pv.sig2_d_init,1));
      end
    end % end getlogPrior
    
  end % end methods 

end % end Dividends