#ifndef __FILE_EMMUSR_H_SEEN__
#define __FILE_EMMUSR_H_SEEN__
/* ----------------------------------------------------------------------------

Copyright (C) 2009.

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 "libsnp.h"
#include "libsmm.h"
#include "emm_base.h"
#include "snp.h"

namespace emm {

  class habit_usrmod;

  typedef habit_usrmod usrmod_type;

  const INTEGER n_state = 1;       //Number of state variables
  const INTEGER n_parms = 8;       //Number of sci mod parameters
  const INTEGER n_stats = 4;       //Number of sci mod functionals
  const INTEGER n_datum = 4;       //Number of simulated variables
  
  class linear_function {
  private:
    REAL a;
    REAL b;
    REAL x0;
  public:
    void initialize(REAL intercept, REAL slope, REAL origin)
      { a = intercept; b = slope; x0 = origin; }
    REAL operator()(REAL x) { return a + b*(x - x0); }
    REAL intercept() { return a; }
    REAL slope() { return b; }
    REAL origin() { return x0; }
  };
  
  class linear_interpolater {
  private:
    std::vector<linear_function> funcs;
    typedef std::vector<linear_function>::size_type fst;
    REAL xmin;
    REAL xmax;
    fst N;
    fst hash(REAL x) { return fst( REAL(N-2)*(x-xmin)/(xmax-xmin) ); }
  public:
    linear_interpolater() 
    {
      scl::realmat grid(2,1) ; grid[1] = 0.0; grid[2] = 1.0;
      scl::realmat vals(2,1) ; vals[1] = 0.0; vals[2] = 1.0;
      update(grid,vals);
    }
    linear_interpolater(scl::realmat& x, scl::realmat& y) { update(x,y); };
    void update(scl::realmat& x, scl::realmat& y) 
    { 
      INTEGER n = x.size(); N = fst(n);  
      funcs.clear(); funcs.reserve(N);
      if (n<2) 
        scl::error("Error, linear_interpolater, x.size() < 2");
      if (x.ncol() != 1 || y.ncol() != 1) 
        scl::error("Error, linear_interpolater, x or y not a vector");
      if (n != y.size()) 
        scl::error("Error, linear_interpolater, x and y sizes differ");
      scl::intvec rank = x.sort(); 
      y = y(rank,"");
      xmin = x[1]; xmax = x[n];
      linear_function f;
      for (INTEGER i=1; i<n; ++i) {
        f.initialize(y[i], (y[i+1]-y[i])/(x[i+1]-x[i]), x[i]);
        funcs.push_back(f);
      }
      funcs.push_back(f);
    }
    REAL operator()(REAL x)
    {
       if (x <= funcs[0].origin()) return funcs[0](x);
       if (x >= funcs[N-1].origin()) return funcs[N-1](x);
       fst i = hash(x);
       if (x < funcs[i].origin()) while(x < funcs[--i].origin());
       else if (x >= funcs[i+1].origin()) while(x >= funcs[++i+1].origin());
       return funcs[i](x);
     }
   };
  
  struct habit_usrmod_variables {
    INTEGER n0;  // for lags and to allow transients to dissipate
    INTEGER n;   // simulation length
    INTEGER n1;  // for leads
    INTEGER nt;  // nt = n0 + n + n1
    scl::realmat log_consumption_growth;
    scl::realmat log_dividend_growth;
    scl::realmat log_surplus_ratio;
    scl::realmat consumption_growth;
    scl::realmat dividend_growth;
    scl::realmat consumption;
    scl::realmat dividend;
    scl::realmat surplus_ratio;
    scl::realmat sensitivity_function;
    scl::realmat marginal_rate_of_substitution;
    scl::realmat price_dividend_ratio;
    scl::realmat bond_price;
    scl::realmat gross_stock_return;
    scl::realmat gross_risk_free_rate;
    scl::realmat geometric_stock_return;
    scl::realmat geometric_risk_free_rate;
    habit_usrmod_variables() { }
    habit_usrmod_variables(INTEGER n_lag, INTEGER n_sim, INTEGER n_lead)
    : n0(n_lag), n(n_sim), n1(n_lead), nt(n0+n+n1),
      log_consumption_growth(nt,1,0.0),
      log_dividend_growth(nt,1,0.0),
      log_surplus_ratio(nt,1,0.0),
      consumption_growth(nt,1,0.0),
      dividend_growth(nt,1,0.0),
      consumption(nt,1,0.0),
      dividend(nt,1,0.0),
      surplus_ratio(nt,1,0.0),
      sensitivity_function(nt,1,0.0),
      marginal_rate_of_substitution(nt,1,0.0),
      price_dividend_ratio(nt,1,0.0),
      bond_price(nt,1,0.0),
      gross_stock_return(nt,1,0.0),
      gross_risk_free_rate(nt,1,0.0),
      geometric_stock_return(nt,1,0.0),
      geometric_risk_free_rate(nt,1,0.0)
    { }
    std::vector<std::string> get_habit_usrmod_variables(scl::realmat& mv);
  };
  
  struct habit_usrmod_parameters {
    REAL g,sig;                // mean and sd. of log consumption growth
    REAL sig_w,rho;            // sd. of log dividend growth and corr. with v
    REAL phi;                  // surplus ratio autoregressive parameter
    REAL delta,gamma;          // time preference and risk aversion parmameters
    scl::realmat R;            // Cholesky factor of variance matix
    REAL mudc;                 // mean of log dividend/consumption series
    INTEGER p;                 // number of model parameters
    habit_usrmod_parameters() 
    : p(n_parms)
    {
      scl::realmat theta(p,1); // These are the values from CC, BGT for mudc.
      theta[1] = 1.57500000000000001e-03; theta[2] = 4.24264068711928486e-03;
      theta[3] = 8.66025403784438695e-04; theta[4] = 3.23316150746190412e-02;
      theta[5] = 9.88461907990555666e-01; theta[6] = 9.90335849608121044e-01;
      theta[7] = 2.00000;                 theta[8] = -3.3965;
      set_theta(theta);
    }
    void set_theta(const scl::realmat& theta)
    {
      if (theta.nrow() != p || theta.ncol() != 1) 
        scl::error("Error, habit_usrmod, set_theta, bad dimension");
      R.resize(2,2);
      g      = theta[1];
      R(1,1) = theta[2];
      R(1,2) = theta[3];
      R(2,2) = theta[4];
      R(2,1) = 0.0;
      scl::realmat sigma = R*T(R);
      sig = sqrt(sigma(1,1));
      sig_w = sqrt(sigma(2,2));
      rho   = sigma(1,2)/(sig*sig_w);
      phi   = theta[5];
      delta = theta[6];
      gamma = theta[7];
      mudc  = theta[8];
    }
    void get_theta(scl::realmat& theta)
    {
      theta.resize(p,1);
      theta[1] = g;
      theta[2] = R(1,1);
      theta[3] = R(1,2);
      theta[4] = R(2,2);
      theta[5] = phi;
      theta[6] = delta;
      theta[7] = gamma;
      theta[8] = mudc;
    }
  };
  
  class habit_usrmod : public libsmm::usrmod_base {
  private:
    habit_usrmod_parameters mp;
    habit_usrmod_variables mv;
    scl::realmat data;
    INT_32BIT bootstrap_seed;
    INT_32BIT simulation_seed;
    linear_interpolater CCpdval;    
    void make_state(INT_32BIT& seed);
    bool make_sim(INT_32BIT& seed, scl::realmat& sim);
  public:
    habit_usrmod(const scl::realmat& dat, INTEGER len_mod_parm,
      INTEGER len_mod_func, const std::vector<std::string>& mod_pfvec,
      const std::vector<std::string>& mod_alvec, std::ostream& detail);
    INTEGER len_rho() { return n_parms; }
    INTEGER len_stats() { return n_stats; }
    bool gen_sim(scl::realmat& sim);
    bool gen_sim(scl::realmat& sim, scl::realmat& func);
    bool gen_bootstrap(std::vector<scl::realmat>& bs);
    void set_rho(const scl::realmat& parm) { mp.set_theta(parm); }
    void get_rho(scl::realmat& parm) { mp.get_theta(parm); }
    bool support(const scl::realmat& theta);
    scl::den_val prior(const scl::realmat& parm, const scl::realmat& func);
    habit_usrmod_variables get_model_variables() { return mv; };
    void write_usrvar(const char* filename);
  };
}

#endif
