/*-----------------------------------------------------------------------------

Copyright (C) 2016

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA   

This program 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 2 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, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-----------------------------------------------------------------------------*/

#include "libscl.h"
#include "crra_mf.h"

using namespace std;
using namespace scl;

// data(1,t) = sr = stock returns as decimal
// data(2,t) = br = bond returns as decimal
// data(3,t) = cg = consumption growth
// data(4,t) = wg = wealth growth (gross return on the consumption endowment)

// Constant Relative Risk Aversion

// mrs = delta*cg^{-gamma}
// lmrs = log(delta) - gamma*lcg

namespace mle {


realmat crra_moment_function::operator()(INTEGER t)
{
  realmat e(2,1);

  REAL lsr = (*dat)(1,t);   //log gross stock return
  REAL lbr = (*dat)(2,t);   //log gross bond return
  REAL lcg = (*dat)(3,t);   //log consumption growth

  REAL lsrlag = (*dat)(1,t-1);
  REAL lbrlag = (*dat)(2,t-1);
  REAL lcglag = (*dat)(3,t-1);

  REAL gamma = mrs_parm[1];
  REAL delta = mrs_parm[2];

  REAL lmrs = log(delta) - gamma*lcg;

  e[1] = 1.0 - exp(lmrs + lsr);
  e[2] = 1.0 - exp(lmrs + lbr);

  realmat m(n_moments,1);

  m[1]  = e[1];
  m[2]  = e[2];
  m[3]  = lsrlag*e[1];
  m[4]  = lbrlag*e[1];
  m[5]  = lcglag*e[1];
  m[6]  = lsrlag*e[2];
  m[7]  = lbrlag*e[2];
  m[8]  = lcglag*e[2];

  return m;
}

realmat crra_moment_function::operator()(INTEGER t, realmat& dm)
{
  realmat e(2,1);

  REAL lsr = (*dat)(1,t);   //log gross stock return
  REAL lbr = (*dat)(2,t);   //log gross bond return
  REAL lcg = (*dat)(3,t);   //log consumption growth

  REAL lsrlag = (*dat)(1,t-1);
  REAL lbrlag = (*dat)(2,t-1);
  REAL lcglag = (*dat)(3,t-1);

  REAL gamma = mrs_parm[1];
  REAL delta = mrs_parm[2];

  REAL lmrs = log(delta) - gamma*lcg;

  e[1] = 1.0 - exp(lmrs + lsr);
  e[2] = 1.0 - exp(lmrs + lbr);

  realmat m(n_moments,1);

  m[1]  = e[1];
  m[2]  = e[2];
  m[3]  = lsrlag*e[1];
  m[4]  = lbrlag*e[1];
  m[5]  = lcglag*e[1];
  m[6]  = lsrlag*e[2];
  m[7]  = lbrlag*e[2];
  m[8]  = lcglag*e[2];

  if (mrs_parm.size() != n_parms) {
    error("Error, crra_moment_function, mrs_parms.size() != n_parms");
  }

  if (n_moments != 8) {
    error("Error, crra_moment_function, n_moments != 8");
  }

  REAL d_lmrs_gamma = -lcg;
  REAL d_lmrs_delta = 1.0/delta;

  REAL d_e1_gamma = - exp(lmrs + lsr)*d_lmrs_gamma;
  REAL d_e2_gamma = - exp(lmrs + lbr)*d_lmrs_gamma;
  REAL d_e1_delta = - exp(lmrs + lsr)*d_lmrs_delta;
  REAL d_e2_delta = - exp(lmrs + lbr)*d_lmrs_delta;

  dm.resize(n_moments,n_parms);

  dm(1,1) = d_e1_gamma;
  dm(2,1) = d_e2_gamma;
  dm(3,1) = lsrlag*d_e1_gamma;
  dm(4,1) = lbrlag*d_e1_gamma;
  dm(5,1) = lcglag*d_e1_gamma;
  dm(6,1) = lsrlag*d_e2_gamma;
  dm(7,1) = lbrlag*d_e2_gamma;
  dm(8,1) = lcglag*d_e2_gamma;

  dm(1,2) = d_e1_delta;
  dm(2,2) = d_e2_delta;
  dm(3,2) = lsrlag*d_e1_delta;
  dm(4,2) = lbrlag*d_e1_delta;
  dm(5,2) = lcglag*d_e1_delta;
  dm(6,2) = lsrlag*d_e2_delta;
  dm(7,2) = lbrlag*d_e2_delta;
  dm(8,2) = lcglag*d_e2_delta;

  return m;
}
}
