#ifndef __FILE_GSMUSR_H_SEEN__
#define __FILE_GSMUSR_H_SEEN__
/* ----------------------------------------------------------------------------

Copyright (C) 2008, 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.

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

//#define DEBUG_GMUSR
#undef DEBUG_GMUSR

/*
Parameter  Name     Meaning                                      Plausible

theta[1]   g        Consumption growth                            0.001575
theta[2]   sig      Standard deviation of log consumption growth  0.00433 
theta[3]   rho      Correlation bewtween lcg and ldg              0.2     
theta[4]   sig_w    Standard deviation of log dividend growth     0.03233 
theta[5]   phi      Autoregressive parameter of surplus_ratio     0.9885 
theta[6]   delta    Subjective discount factor                    0.9903 
theta[7]   gamma    Risk aversion parameter                       2.00

stats[1]   br       Bond returns                                  1.000743618
stats[2]   bv       Variance of bond returns                      0.000024083
stats[3]   sr       Stock returns                                 1.002851
stats[4]   sv       Variance of stock returns                     0.001771989

sig^2   = R11^2 + R12^2                                           0.00433^2
rho     = R12*R22/(sig*sig_w)                                     0.2
sig_w^2 = R22^2                                                   0.03233^2

R11 = sig*sqrt(1 - rho^2)                                         0.004243
R12 = rho*sig                                                     0.000866
R22 = sig_w                                                       0.032332

Parameter  Name     Reasonable support

theta[1]   g        [-0.02, 0.02]
theta[2]   sig      [ 0.000, 1.000]
theta[3]   rho      [-1.000, 1.000]
theta[4]   sig_w    [ 0.000, 3.000]
theta[5]   phi      [-1.000, 1.000]              
theta[6]   delta    [0.90, 0.99999]
theta[7]   gamma    [0.25, 10.0]
*/

#include "libgsm.h"
#include "libsnp.h"
#include "snp.h"

namespace gsm {

  class habit_sci_mod;
  class stat_mod;
  class snp_stat_mod;

  typedef habit_sci_mod sci_mod_type;
  typedef snp_stat_mod stat_mod_type;

  const INTEGER n_state = 1;       //Number of state variables
  const INTEGER n_sci_parms = 7;   //Number of sci mod parameters
  const INTEGER n_sci_funcs = 4;   //Number of sci mod functionals
  const INTEGER n_stat_funcs = 9;  //Number of stat mod functionals
  
  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 lfst;
    REAL xmin;
    REAL xmax;
    lfst N;
    lfst hash(REAL x) { return lfst( 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 = lfst(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);
       lfst 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_sci_mod_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 dividend;
    scl::realmat consumption;
    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_sci_mod_variables() { }
    habit_sci_mod_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),
      dividend(nt,1,0.0),
      consumption(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_sci_mod_variables(scl::realmat& mv);
  };
  
  struct habit_sci_mod_parameters {
    REAL g;                // Consumption growth
    REAL sig;              // Standard deviation of log consumption growth
    REAL rho;              // Correlation bewtween lcg and ldg 
    REAL sig_w;            // Standard deviation of log dividend growth
    REAL phi;              // Autoregressive parameter of surplus_ratio
    REAL delta;            // Subjective discount factor
    REAL gamma;            // Risk aversion parameter
    scl::realmat R;        // Cholesky factor of variance matix
    const INTEGER p;       // Number of model parameters
    habit_sci_mod_parameters() 
    : p(n_sci_parms)
    {
      scl::realmat theta(p,1); // These are the values from CC.
      theta[1] = 1.57500000000000001e-03; 
      theta[2] = 4.33012701892219330e-03;
      theta[3] = 2.00000000000000000e-01; 
      theta[4] = 3.23316150746190412e-02;
      theta[5] = 9.88461907990555666e-01; 
      theta[6] = 9.90335849608121044e-01;
      theta[7] = 2.00000;
      set_theta(theta);
    }
    void set_theta(const scl::realmat& theta)
    {
      if (theta.nrow() != p || theta.ncol() != 1) 
        scl::error("Error, habit_sci_mod, set_theta, bad dimension");
      R.resize(2,2);
      g      = theta[1];
      sig    = theta[2];
      rho    = theta[3];
      sig_w  = theta[4];
      phi    = theta[5];
      delta  = theta[6];
      gamma  = theta[7];
      R(1,1) = sig*sqrt(1.0 - pow(rho,2));
      R(2,1) = 0.0;
      R(1,2) = rho*sig;
      R(2,2) = sig_w;
    }
    void get_theta(scl::realmat& theta)
    {
      theta.resize(p,1);
      theta[1] = g;
      theta[2] = sig;
      theta[3] = rho;
      theta[4] = sig_w;
      theta[5] = phi;
      theta[6] = delta;
      theta[7] = gamma;
    }
  };

  class habit_sci_mod : public libgsm::sci_mod_base {
  private:
    habit_sci_mod_parameters mp;
    habit_sci_mod_variables mv;
    void make_state(INT_32BIT& seed);
    bool gen_sim_cgsrpd(scl::realmat& sim);
    bool gen_sim_cgsrpd(scl::realmat& sim, scl::realmat& func);
    linear_interpolater CCpdval;    
    std::ostream& debug;
    bool debug_switch;
    INTEGER n_dat_vars;    //Number of simulated variables
  public:
    habit_sci_mod
      (const scl::realmat* dat_ptr, const std::vector<std::string>& pfvec,
       const std::vector<std::string>& alvec, std::ostream& detail);
    INTEGER len_parm() { return n_sci_parms; }
    INTEGER len_func() { return n_sci_funcs; }
    void get_parm(scl::realmat& parm) { mp.get_theta(parm); }
    void set_parm(const scl::realmat& parm) { mp.set_theta(parm); }
    bool support(const scl::realmat& parm);
    bool gen_sim(scl::realmat& sim, scl::realmat& func);
    scl::den_val prior(const scl::realmat& parm, const scl::realmat& func);
    scl::realmat annualize_parms(const scl::realmat& theta);
    scl::realmat annualize_funcs(const scl::realmat& stats);
  };
  
  class stat_mod : public libgsm::stat_mod_base {
  private:
    scl::realmat  b, B; // location function parameters
    scl::realmat  p, P; // variance function parameters, the ARCH or MA part
    scl::realmat  G;    // variance function parameters, the GARCH or AR part
    const scl::realmat* Y;    // the data, which is Y(d,n)
    INTEGER d;  // dimension d of the model, inferred from row dimension of Y
    INTEGER n;  // sample size n, inferred from column dimension of Y
    INTEGER lagu; // lags in mean function
    INTEGER lagR; // lags in MA-part of variance function
    INTEGER lagG; // lags in AR-part of variance function
    INTEGER lR; // dimension of vech(R), R is Cholesky factor of variance
    INTEGER len_sim;
    REAL get_eta(INTEGER i);
    void set_eta(INTEGER i, REAL v);
  public:
    stat_mod
      (const scl::realmat* dat_ptr, const std::vector<std::string>& pfvec,
       const std::vector<std::string>& alvec, std::ostream& detail);
    INTEGER len_parm() { return d + d*d*lagu + lR + lR*d*lagR + lagG; } 
    INTEGER len_func() { return n_stat_funcs; }
    void get_parm(scl::realmat& parm);
    void set_parm(const scl::realmat& parm);
    void set_data_ptr(const scl::realmat* dat_ptr);
    bool gen_sim(scl::realmat& sim, scl::realmat& func);
    bool support(const scl::realmat& parm);
    scl::den_val prior(const scl::realmat& parm, const scl::realmat& func);
    scl::den_val logl();
  };

  class snp_stat_mod : public libgsm::stat_mod_base {
  private:
    scl::realmat Xlags;      // First few transformed and squashed
    scl::realmat Ylags;      // observations to get lags for gen_sim
    scl::realmat X;          // Transformed and squashed data or simulation
    scl::realmat Y;          // maintained by set_data_ptr for use by logl, 
    snp::trnfrm tr;
    snp::snpll ll;
    INTEGER lsim;
    INTEGER spin;
    INTEGER lrho;
    INTEGER lstats;
    bool snp_simulate(INT_32BIT& iseed, scl::realmat& sim);
  public:
    snp_stat_mod
      (const scl::realmat* dat_ptr, const std::vector<std::string>& pfvec,
       const std::vector<std::string>& alvec, std::ostream& detail);
    INTEGER len_parm() {return lrho;}
    INTEGER len_func() {return lstats;}
    INTEGER num_obs() {return Y.get_cols();}
    void    get_parm(scl::realmat& rho);
    void    set_parm(const scl::realmat& rho);
    void    set_data_ptr(const scl::realmat* dat_ptr);
    bool    support(const scl::realmat& rho);
    scl::den_val prior(const scl::realmat& rho, const scl::realmat& stats);
    bool    gen_sim(scl::realmat& sim, scl::realmat& stats);
    scl::den_val logl();
    scl::den_val logl(scl::realmat& dlogl);
    REAL    penalty(scl::realmat& dpenalty);
  };
}

#endif
