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

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.

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

#include "gsmusr.h"
#include "libsnp.h"

using namespace std;
using namespace scl;

namespace gsm {

  vector<string> 
  habit_sci_mod_variables::get_habit_sci_mod_variables(realmat& mv)
  {
    vector<string> names;
    names.push_back("log_consumption_growth");
    mv = log_consumption_growth;
    names.push_back("log_dividend_growth");
    mv = cbind(mv,log_dividend_growth);
    names.push_back("log_surplus_ratio");
    mv = cbind(mv,log_surplus_ratio);
    names.push_back("consumption_growth");
    mv = cbind(mv,consumption_growth);
    names.push_back("dividend_growth");
    mv = cbind(mv,dividend_growth);
    names.push_back("consumption");
    mv = cbind(mv,consumption);
    names.push_back("dividend");
    mv = cbind(mv,dividend);
    names.push_back("surplus_ratio");
    mv = cbind(mv,surplus_ratio);
    names.push_back("sensitivity_function");
    mv = cbind(mv,sensitivity_function);
    names.push_back("marginal_rate_of_substitution");
    mv = cbind(mv,marginal_rate_of_substitution);
    names.push_back("price_dividend_ratio");
    mv = cbind(mv,price_dividend_ratio);
    names.push_back("bond_price"); 
    mv = cbind(mv,bond_price);
    names.push_back("gross_stock_return");
    mv=cbind(mv,gross_stock_return);
    names.push_back("gross_risk_free_rate");
    mv=cbind(mv,gross_risk_free_rate);
    names.push_back("geometric_stock_return");
    mv=cbind(mv,geometric_stock_return);
    names.push_back("geometric_risk_free_rate");
    mv=cbind(mv,geometric_risk_free_rate);
    realmat log_price_dividend_ratio = price_dividend_ratio;
    for (INTEGER i=1; i<=3; ++i) log_price_dividend_ratio[i] = 0.0;
    for (INTEGER i=4; i<=nt; ++i) 
      log_price_dividend_ratio[i] = log(log_price_dividend_ratio[i]);
    names.push_back("log_price_dividend_ratio");
    mv = cbind(mv,log_price_dividend_ratio);
    return names;
  }

  habit_sci_mod::habit_sci_mod
    (const scl::realmat* dat_ptr, const std::vector<std::string>& pfvec,
     const std::vector<std::string>& alvec, std::ostream& detail)
  :  debug(detail), debug_switch(false)
  {
    n_dat_vars = dat_ptr->nrow();

    vector<string>::const_iterator alvec_ptr = alvec.begin();

    INTEGER n0 = atoi((++alvec_ptr)->substr(0,12).c_str());
    INTEGER n  = atoi((++alvec_ptr)->substr(0,12).c_str());
    INTEGER n1 = atoi((++alvec_ptr)->substr(0,12).c_str());

    mv = habit_sci_mod_variables(n0,n,n1);

    realmat sgrid(22,1), pdval(22,1);  // CC policy function
    sgrid[1] = -5.20514;  pdval[1] = 4.35107;
    sgrid[2] = -4.96183;  pdval[2] = 4.44397;
    sgrid[3] = -4.76405;  pdval[3] = 4.52012;
    sgrid[4] = -4.74214;  pdval[4] = 4.5287;
    sgrid[5] = -4.52244;  pdval[5] = 4.61509;
    sgrid[6] = -4.27914;  pdval[6] = 4.7138;
    sgrid[7] = -4.07091;  pdval[7] = 4.80184;
    sgrid[8] = -3.66544;  pdval[8] = 4.9827;
    sgrid[9] = -3.37776;  pdval[9] = 5.11901;
    sgrid[10] = -3.15462;  pdval[10] = 5.23134;
    sgrid[11] = -2.97229;  pdval[11] = 5.32856;
    sgrid[12] = -2.86453;  pdval[12] = 5.38912;
    sgrid[13] = -2.81814;  pdval[13] = 5.41611;
    sgrid[14] = -2.68461;  pdval[14] = 5.49662;
    sgrid[15] = -2.56683;  pdval[15] = 5.57113;
    sgrid[16] = -2.46147;  pdval[16] = 5.64058;
    sgrid[17] = -2.40616;  pdval[17] = 5.67804;
    sgrid[18] = -2.39616;  pdval[18] = 5.6849;
    sgrid[19] = -2.38616;  pdval[19] = 5.69179;
    sgrid[20] = -2.37616;  pdval[20] = 5.69873;
    sgrid[21] = -2.37191;  pdval[21] = 5.7017;
    sgrid[22] = -2.36616;  pdval[22] = 5.70575;

    CCpdval.update(sgrid,pdval);
    
    #if defined DEBUG_GMUSR

      debug_switch = true;

      detail << starbox("/values read from parmfile//");
      detail << '\n';
      detail << "\t\t\t\t n0 = " << n0 << '\n';
      detail << "\t\t\t\t n = " << n << '\n';
      detail << "\t\t\t\t n1 = " << n1 << '\n';
      detail.flush();

      detail << starbox("/CC monthly variance matrix and its factor//");
      realmat V(2,2);
      V(1,1) = 1.5/(100.0*sqrt(12));
      V(2,2) = 11.2/(100.0*sqrt(12));
      V(1,2) = 0.2*V(1,1)*V(2,2);
      V(2,1) = V(1,2);
      V(1,1) = pow(V(1,1),2);
      V(2,2) = pow(V(2,2),2);
      realmat A(3,1), J(3,3);
      A[1] = V(1,1);
      A[2] = V(1,2);
      A[3] = V(2,2);
      libsnp::factor(A.get_x(),J.get_x(),2);
      detail << V << A;
      detail << starbox("/theta values from CC//") << '\n';
      detail << '\n';
      detail << "\t\t\t" << fmt('e',27,17,1.89/(100.0*12)) << '\n';
      detail << "\t\t\t" << fmt('e',27,17,1.5/(100.0*sqrt(12.0))) << '\n';
      detail << "\t\t\t" << fmt('e',27,17,0.2) << '\n';
      detail << "\t\t\t" << fmt('e',27,17,11.2/(100.0*sqrt(12.0))) << '\n';
      detail << "\t\t\t" << fmt('e',27,17,pow(0.87,1.0/12.0)) << '\n';
      detail << "\t\t\t" << fmt('e',27,17,pow(0.89,1.0/12.0)) << '\n';
      detail << "\t\t\t" << fmt('e',27,17,2.0) << '\n';
      
      realmat theta;
      mp.get_theta(theta);
      detail << starbox("/theta from mp struct//");
      detail << theta;
      detail << starbox("/implied monthly and annual CC values//");
      detail << '\n';
      detail << "\t\t\t\t" << "g       " << mp.g     << '\n';
      detail << "\t\t\t\t" << "sig     " << mp.sig   << '\n';
      detail << "\t\t\t\t" << "sig_w   " << mp.sig_w << '\n';
      detail << "\t\t\t\t" << "rho     " << mp.rho   << '\n';
      detail << "\t\t\t\t" << "phi     " << mp.phi   << '\n';
      detail << "\t\t\t\t" << "delta   " << mp.delta << '\n';
      detail << "\t\t\t\t" << "gamma   " << mp.gamma << '\n';
      detail << '\n';
      detail << "\t\t\t\t" << "g       " << 100.0*12.0*mp.g << '\n';
      detail << "\t\t\t\t" << "sig     " << 100.0*sqrt(12)*mp.sig << '\n';
      detail << "\t\t\t\t" << "sig_w   " << 100.0*sqrt(12)*mp.sig_w << '\n';
      detail << "\t\t\t\t" << "rho     " << mp.rho << '\n';
      detail << "\t\t\t\t" << "phi     " << pow(mp.phi,12) << '\n';
      detail << "\t\t\t\t" << "delta   " << pow(mp.delta,12) << '\n';
      detail << "\t\t\t\t" << "gamma   " << mp.gamma << '\n';
      detail.flush();

      realmat sim, stats;
      gen_sim_cgsrpd(sim,stats);

      realmat variables;
      vector<string> names = mv.get_habit_sci_mod_variables(variables);

      INTEGER r = variables.nrow();
      INTEGER c = variables.ncol();

      realmat sum(c,1);
      
      fill(sum,0.0);
      for (INTEGER j=1; j<=c; ++j) {
        for (INTEGER i=mv.n0+1; i<=mv.n0+mv.n; ++i) {
          sum[j] += variables(i,j);
        }
      }
      realmat mean = sum/mv.n;

      fill(sum,0.0);
      for (INTEGER j=1; j<=c; ++j) {
        for (INTEGER i=mv.n0+1; i<=mv.n0+mv.n; ++i) {
          sum[j] += pow(variables(i,j)-mean[j],2);
        }
      }
      realmat var = sum/mv.n;

      detail << starbox("/means and standard errors of model variables//");
      detail << '\n';
      for (INTEGER j=1; j<=c; ++j) {
        if (names[j-1] == "consumption") {
          detail << '\t' << names[j-1] << ' ' 
            << fmt('e',45-names[j-1].size(),5,mean[j]) 
            << fmt('e',15,5,sqrt(var[j])) << '\n';
         }
         else {
          detail << '\t' << names[j-1] << ' ' 
            << fmt('f',45-names[j-1].size(),5,mean[j]) 
            << fmt('f',15,5,sqrt(var[j])) << '\n';
         }
      }
      detail.flush();

      r = sim.nrow();
      c = sim.ncol();
      sum.resize(r,1);

      fill(sum,0.0);
      for (INTEGER j=1; j<=c; ++j) {
        for (INTEGER i=1; i<=r; ++i) {
          sum[i] += sim(i,j);
        }
      }
      mean = sum/c;

      fill(sum,0.0);
      for (INTEGER j=1; j<=c; ++j) {
        for (INTEGER i=1; i<=r; ++i) {
          sum[i] += pow(sim(i,j)-mean[i],2);
        }
      }
      var = sum/c;

      detail << starbox("/means and standard errors of annual simulation//");
      detail << '\n';
      detail << "\t\t" << "consumption growth  " 
             << fmt('f',9,5,mean[1]) <<  fmt('f',9,5,sqrt(var[1])) << '\n'; 
      detail << "\t\t" << "stock returns       " 
             << fmt('f',9,5,mean[2]) <<  fmt('f',9,5,sqrt(var[2])) << '\n'; 

      debug_switch = false;

    #endif
    
  }
  
  void habit_sci_mod::make_state(INT_32BIT& seed)
  {
    const INT_32BIT start_seed = 100542;
  
    REAL del_c;
    REAL del_d;
    
    REAL cg;
    REAL dg;
    
    REAL s; 
    REAL S_bar = mp.sig*sqrt(mp.gamma/(1.0-mp.phi));
    REAL s_bar = log(S_bar);
    REAL s_max = s_bar + (0.5)*(1.0 - pow(S_bar,2));
    REAL s_lag = s_bar;
    REAL del_s;
  
    REAL lam; 
    REAL lam_lag = 0.0;
  
    REAL mrs;
    
    realmat e(2,1);
    
    INT_32BIT state_seed = start_seed;
  
    REAL C = REAL_MIN;
    REAL D = REAL_MIN;

    for (INTEGER t=2; t <=mv.nt; t++) {
  
      if (t == mv.n0+1) state_seed = seed;
  
      e[1] = unsk(state_seed);
      e[2] = unsk(state_seed);
  
      e = mp.R*e;
   
      del_c = mp.g + e[1];
      del_d = mp.g + e[2];
  
      cg = exp(del_c);
      dg = exp(del_d);
    
      C *= cg;
      D *= dg;
  
      mv.log_consumption_growth[t] = del_c;
      mv.consumption_growth[t] = cg;
      mv.consumption[t] = C;
      
      mv.log_dividend_growth[t] = del_d;
      mv.dividend_growth[t] = dg;
      mv.dividend[t] = D;
      
      s = (1.0 - mp.phi)*s_bar + mp.phi*s_lag + lam_lag*e[1];
      del_s = s - s_lag;
      s_lag = s;
      
      mv.log_surplus_ratio[t] = s;
      mv.surplus_ratio[t] = exp(s);
  
      lam = (s < s_max) ? sqrt(1.0 - 2.0*(s-s_bar))/S_bar - 1.0 : 0.0;
      lam_lag = lam;
      
      mv.sensitivity_function[t] = lam;
  
      mrs = mp.delta*exp(-mp.gamma*(del_s + del_c));
  
      mv.marginal_rate_of_substitution[t] = mrs;
    }
    seed = state_seed;
    if (!IsFinite(C)) error("Error, sci_mod::make_state, C not finite");
  }
  
  bool habit_sci_mod::gen_sim_cgsrpd(realmat& sim)
  {
    realmat theta(mp.p,1);
    this -> get_parm(theta);
      
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << starbox("/Entered gen_sim_cgsrpd with these parameter values//");
      debug << theta << '\n';
    }
    #endif

    INT_32BIT sim_seed = 1237173085;
    this->make_state(sim_seed);

    REAL smax = -REAL_MAX;
    for (INTEGER t=mv.n0+1; t<=mv.n0+mv.n; ++t) {
      smax = smax < mv.log_surplus_ratio[t] ? mv.log_surplus_ratio[t] : smax;
    }

    REAL S_bar = mp.sig*sqrt(mp.gamma/(1.0-mp.phi));
    REAL s_bar = log(S_bar);
    REAL s_max = s_bar + (0.5)*(1.0 - pow(S_bar,2));
    REAL S_max = exp(s_max);

    smax = smax < s_max ? s_max : smax;

    REAL top = exp(smax);

    const INTEGER intervals = 10;        // the following construction of
                                         // sgrid replicates that in CC
    REAL increment = top/(intervals+1);
    
    bool no_S_bar = true;
    bool no_S_max = true;;
    realmat sgrid;
    for (INTEGER i=1; i<=intervals; ++i) {  // equal increments in S
      REAL Si = i*increment;
      if (Si == S_bar) no_S_bar = false;
      if (Si == S_max) no_S_max = false;
      sgrid.push_back(log(Si));
    }
    if (no_S_bar) sgrid.push_back(s_bar);  // make sure these two in grid
    if (no_S_max) sgrid.push_back(s_max);
    sgrid.push_back(smax - 0.01);          // extra points added near smax
    sgrid.push_back(smax - 0.02);
    sgrid.push_back(smax - 0.03);
    sgrid.push_back(smax - 0.04);

    sgrid.sort();
    REAL smin = sgrid[1];

    const INTEGER npoints = 5;    // adding quadrature points to sgrid
                                  // is a departure from CC
    realmat x, w;
    if( hquad(npoints,x,w) ) error("Error, sci_mod, gen_sim_cgsrpd, hquad failed");
    
    const REAL root2 = sqrt(2.0);              // 1.4142135623730951
    const REAL rootpi = sqrt(4.0*atan(1.0));   // 1.7724538509055161
    
    realmat abscissae(npoints,1);
    realmat weights(npoints,1);
    for (INTEGER i=1; i<=npoints; ++i) {
      abscissae[i] = root2*mp.sig*x[i];
      weights[i] = w[i]/rootpi;
    }

    for (INTEGER j=1; j<=npoints; ++j) {
      REAL s0 = smin;
      REAL lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
      REAL e1 = abscissae[j];
      REAL s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
      sgrid.push_back(s1);
      s0 = smax;
      lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
      e1 = abscissae[j];
      s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
      sgrid.push_back(s1);
    }
    
    sgrid.sort();

    INTEGER ngrid = sgrid.size();          // delete points that are too close
    intvec idx = seq(1,ngrid);
    for (INTEGER i=1; i<ngrid; ++i) {
      if (sgrid[i+1] - sgrid[i] < 0.001) idx[i+1] = -1;
    }

    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << starbox("/Check interpolater//") << '\n';
      debug << '\n';
      debug << "\t\t" << "ngrid before purge = " << ngrid << '\n';
    }
    #endif

    sgrid = sgrid(idx,"");
    ngrid = sgrid.size();

    intvec tstidx = seq(1,ngrid);       // only check convergence between
    INTEGER mid = (ngrid + 1)/2;        // smin and smax
    for (INTEGER i=1; i<=ngrid; ++i) {
      if (sgrid[i] < smin || smax < sgrid[i]) tstidx[i] = mid;
    }

    realmat pdval(ngrid,1);

    for (INTEGER i=1; i<=ngrid; ++i) {
      pdval[i] = CCpdval(sgrid[i]);
    }
    
    linear_interpolater v(sgrid,pdval);  
    
    // v maps log_surplus_ratio to log_price_dividend_ratio, i.e.,
    // log_price_dividend_ratio = v(log_surplus_ratio)
    
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << "\t\t" << "ngrid after purge = " << ngrid << '\n';
      debug << '\n';
      debug << "\t\t" << "sgrid";
      debug << '\n';
      debug << sgrid;
      debug << '\n';
      debug << "\t\t" << "pdval";
      debug << pdval;
      for (INTEGER i=1; i<=ngrid; ++i) { 
        pdval[i] =  36 + 5*sgrid[i] + exp(sgrid[i]);
      }
      debug << '\n';
      debug << "\t\t" << "modified pdval";
      debug << pdval;
      debug << '\n';
      debug << "\t\t" << "s, pd_true, pd_approx, error" << '\n';
      debug << '\n';
      v.update(sgrid,pdval);
      for (INTEGER i=2; i<=101; i+=5) {
        REAL x = log(i*S_max/101.0);
        REAL y1 =  36 + 5*x+ exp(x);
        REAL y2 = v(x);
        debug << "\t\t" <<fmt('f',10,5,x) << fmt('f',10,5,y1) 
          << fmt('f',10,5,y2) << fmt('f',10,5,100.0*(y2-y1)/y1) <<'\n';
      }
      for (INTEGER i=1; i<=ngrid; ++i) {
        pdval[i] = CCpdval(sgrid[i]);
      }
      v.update(sgrid,pdval);
      debug << '\n';
      debug << "\t\t" << "reset pdval";
      debug << pdval;
    }
    #endif 

    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << starbox("/quadrature points//");
      debug << cbind(abscissae,weights);
    }
    #endif
    
    const INTEGER maxiter = 2000;
    const REAL tol = 1.0e-4;
    bool converge = false;

    realmat PD0(ngrid,1);
    realmat PD0_lag(ngrid,1);

    for (INTEGER i=1; i<=ngrid; ++i) {
      PD0_lag[i] = exp(v(sgrid[i]));
    }

    #if defined DEBUG_GMUSR
      string pfilename = "hab.policy_function_vals.dat";
      string qfilename = "hab.policy_function_qpts.dat";
      string sfilename = "hab.simulation.dat";
      realmat policy = sgrid;
      policy = cbind(policy,pdval);
      realmat qpts(2*npoints,1);
      realmat qval(2*npoints,1);
      for (INTEGER j=1; j<=npoints; ++j) {
        REAL s0 = smin;
        REAL lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
        REAL e1 = abscissae[j];
        REAL s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
        qpts[j] = s1;
        qval[j] = v(s1);
        s0 = smax;
        lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
        e1 = abscissae[j];
        s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
        qpts[npoints+j] = s1;
        qval[npoints+j] = v(s1);
      }
      qpts = cbind(qpts,qval);
    #endif
    
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      string msg = "/Policy function iterations written to/";
      msg += pfilename + "/" + qfilename + "//";
      msg += "Annual simulation written to/";
      msg += sfilename + "//";
      debug << starbox(msg.c_str());
    }
    #endif

    REAL Ewv_fixed_part = exp(0.5*(1.0 - pow(mp.rho,2))*pow(mp.sig_w,2));
    REAL G = exp(mp.g);

    for (INTEGER iter=1; iter<=maxiter; ++iter) {
      for (INTEGER i=1; i<=ngrid; ++i) {
        REAL s0 = sgrid[i];
        REAL lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
        PD0[i] = 0.0;
        for (INTEGER j=1; j<=npoints; ++j) {
          REAL e1 = abscissae[j];
          REAL s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
          REAL del_s = s1 - s0;
          REAL del_c = mp.g + e1;
          REAL M1 = mp.delta*exp(-mp.gamma*(del_s + del_c));
          REAL PD1 = exp(v(s1));
          REAL Ewv = Ewv_fixed_part*exp(mp.rho*(mp.sig_w/mp.sig)*e1);
          PD0[i] += G*M1*Ewv*(1.0 + PD1)*weights[j];
        }
      }
  
      for (INTEGER i=1; i<=ngrid; ++i) {
        pdval[i] = log(PD0[i]);
      }
      v.update(sgrid,pdval);

      #if defined DEBUG_GMUSR
      if (debug_switch) {
        policy = cbind(policy,pdval);
        for (INTEGER i=1; i<=2*npoints; ++i) {
          REAL s1 = qpts(i,1);
          qval[i] = v(s1);
        }
        qpts = cbind(qpts,qval);
      }
      #endif

      converge = true;
      for (INTEGER i=1; i<=ngrid; ++i) {
        INTEGER j = tstidx[i];
        if (fabs(PD0[j] - PD0_lag[j]) > fabs(PD0_lag[j] + tol)*tol) {
          converge = false;
          break;
        }
      }

      #if defined DEBUG_GMUSR
      if (debug_switch) {
        debug << "\t\t" << "iter = " << iter << ' ';
        debug << "converge = " << boolalpha <<' '<< converge <<'\n';
      }
      #endif

      for (INTEGER i=1; i<=ngrid; ++i) {
        PD0_lag[i] = PD0[i];
      }

      if (converge) break;
    }

    if (!converge) {
      warn("Warning, convergence failed at this parameter value");
      cerr << theta;
    }

    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << starbox("/log price dividend ratio interpolater//");
      debug << '\n';;
      for (INTEGER i=1; i<=ngrid; ++i) {
      debug << '\t' << "sgrid[" << i << "] = " << sgrid[i] << ";  "
           << "pdval[" << i << "] = " << pdval[i] << ";\n";
      }
      converge = false;
    }
    #endif
    
    realmat bval(ngrid,1,-0.00119638);    // fill correspondes to 
                                          // CC parameter values
    linear_interpolater b(sgrid,bval);
    
    // b maps log_surplus_ratio to the log_bond_price, i.e.,
    // log_bond_price = B(log_surplus_ratio)

    realmat B0(ngrid,1);
    
    for (INTEGER i=1; i<=ngrid; ++i) {
      REAL s0 = sgrid[i];
      REAL lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
      B0[i] = 0.0;
      for (INTEGER j=1; j<=npoints; ++j) {
        REAL e1 = abscissae[j];
        REAL s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
        REAL del_s = s1 - s0;
        REAL del_c = mp.g + e1;
        REAL M1 = mp.delta*exp(-mp.gamma*(del_s + del_c));
        B0[i] += M1*weights[j];
      }
    }

    for (INTEGER i=1; i<=ngrid; ++i) {
      bval[i] = log(B0[i]);
    }
    
    b.update(sgrid,bval);

    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << starbox("/bond log price interpolater/sgrid, bval//");
      debug << cbind(sgrid,bval);
    }
    #endif

    for (INTEGER t=2; t <= mv.nt; ++t) {
      mv.price_dividend_ratio[t] = exp(v(mv.log_surplus_ratio[t]));
      mv.bond_price[t] = exp(b(mv.log_surplus_ratio[t]));
    }

    for (INTEGER t=3; t <= mv.nt; ++t) {
      REAL ret = 1.0 + mv.price_dividend_ratio[t];
      ret /= mv.price_dividend_ratio[t-1];
      ret *= mv.dividend_growth[t];
      mv.gross_stock_return[t] = ret;
      mv.gross_risk_free_rate[t] = (1.0)/mv.bond_price[t];
      mv.geometric_stock_return[t] = log(ret);
      mv.geometric_risk_free_rate[t] = -log(mv.bond_price[t]);
    }

    sim.resize(3,mv.n/12);

    REAL r0,C0,C12,D0,D12,PD; 
    INTEGER tt = 0;
    for (INTEGER t=1; t <= mv.n/12; ++t) {
      r0 = C0 = C12 = D0 = D12 = 0.0;
      for (INTEGER j=1; j <= 12; j++) {
        ++tt;
        r0  += mv.geometric_stock_return[mv.n0 + tt];
        C0  += mv.consumption[mv.n0 + tt];
        C12 += mv.consumption[mv.n0 + tt - 12];
        D0  += mv.dividend[mv.n0 + tt];
        D12 += mv.dividend[mv.n0 + tt - 12];
      }
      PD = mv.price_dividend_ratio[mv.n0 + tt]*mv.dividend[mv.n0 + tt]/D0;
      sim(1,t) = log(C0) - log(C12);
      sim(2,t) = r0;
      sim(3,t) = log(PD);
    }

    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << starbox("/First 24 simulated monthly values//");
      realmat monthly;
      vector<string> names = mv.get_habit_sci_mod_variables(monthly);
      idx = seq(mv.n0+1, mv.n0+24);
      monthly = monthly(idx,"");
      for (vector<string>::size_type i=0; i<names.size(); ++i ){ 
        debug << '\n';
        debug << names[i] << monthly("",i+1);;
      }
      debug << starbox("/First 5 simulated annual values//");
      idx = seq(1,5);
      debug << sim("",idx);

      vecwrite(pfilename.c_str(), policy);
      vecwrite(qfilename.c_str(), qpts);
      vecwrite(sfilename.c_str(), sim);
    }
    #endif    
    
    return converge; 
  }
  
  bool habit_sci_mod::gen_sim_cgsrpd(realmat& sim, realmat& stats)
  {
    bool success = this->gen_sim_cgsrpd(sim);
  
    if (success) {
    
      REAL mean_rf = 0.0;
      REAL mean_sr = 0.0;
  
      for (INTEGER t=mv.n0+1; t <= mv.n0+mv.n; t++) {
        mean_rf += mv.geometric_risk_free_rate[t];
        mean_sr += mv.geometric_stock_return[t];
      }
  
      mean_rf /= REAL(mv.n);
      mean_sr /= REAL(mv.n);
    
      REAL var_rf = 0.0;
      REAL var_sr = 0.0;
      
      for (INTEGER t=mv.n0+1; t <= mv.n0+mv.n; t++) {
        var_rf += pow(mv.geometric_risk_free_rate[t]-mean_rf, 2);
        var_sr += pow(mv.geometric_stock_return[t]-mean_sr, 2);
      }
  
      var_rf /= REAL(mv.n);
      var_sr /= REAL(mv.n);
    
      stats.resize(n_sci_funcs,1);
  
      #if defined DEBUG_GMUSR
        stats.check(1) = mean_rf;
        stats.check(2) = var_rf;
        stats.check(3) = mean_sr;
        stats.check(4) = var_sr;
      #else
        stats[1] = mean_rf;
        stats[2] = var_rf;
        stats[3] = mean_sr;
        stats[4] = var_sr;
      #endif

    }
    return success;
  }

  bool habit_sci_mod::support(const realmat& theta) 
  {

    realmat theta_lo(n_sci_parms,1);
    
    theta_lo[1] = -0.02;      // g 
    theta_lo[2] =  0.0;       // sig   
    theta_lo[3] = -1.0;       // rho   
    theta_lo[4] =  0.0;       // sig_w 
    theta_lo[5] = -1.0;       // phi  
    theta_lo[6] =  0.9;       // delta
    theta_lo[7] =  0.25;      // gamma
  
    realmat theta_hi(n_sci_parms,1);
  
    theta_hi[1] =  0.02;      // g 
    theta_hi[2] =  1.0;       // sig   
    theta_hi[3] =  1.0;       // rho   
    theta_hi[4] =  3.0;       // sig_w 
    theta_hi[5] =  1.0;       // phi  
    theta_hi[6] =  0.99999;   // delta
    theta_hi[7] = 10.0;       // gamma
  
    return ( (theta_lo < theta) && (theta < theta_hi) );
  }

  den_val habit_sci_mod::prior(const realmat& theta,const realmat& stats)
  {
    const REAL minus_log_root_two_pi = -9.1893853320467278e-01;
  
    const REAL bond = 0.000743618;             // 0.89% annualized
    const REAL sbond = (1.0/1200.0)/1.96;      // P(within 1% annualized)=0.95
  
    realmat mean(n_sci_parms,1,0.0);
    realmat sdev(n_sci_parms,1,100.0);

    mean[ 1] =   0.001575;      // g
    sdev[ 1] =   0.001/1.96;
    mean[ 2] =   0.00433013;    // sig
    sdev[ 2] =   0.0001/1.96;
    mean[ 3] =   0.2;           // rho 
    sdev[ 3] =   0.1/1.96;
    mean[ 4] =   0.0323316;     // sig_w
    sdev[ 4] =   0.001/1.96;
    mean[ 5] =   0.988462;      // phi
    sdev[ 5] =   0.01/1.96;
    mean[ 6] =   0.990336;      // delta
    sdev[ 6] =   0.001/1.96;
    mean[ 7] =   2.0;           // gamma
    sdev[ 7] =   1.0/1.96;

    for (INTEGER i=1; i<=n_sci_parms; ++i) {
      sdev[i] = fabs(mean[i])*0.1/1.96;
    }

    sdev[5] = fabs(mean[5])*0.001/1.96;    // Don't forget to report
    sdev[6] = fabs(mean[6])*0.001/1.96;    // this in the paper.

    den_val sum(true,0.0);
    
    REAL zbond = (stats[1]-bond)/sbond;
    REAL ebond = minus_log_root_two_pi - log(sbond) - 0.5*pow(zbond,2);
  
    sum += den_val(true,ebond);
  
    for (INTEGER i=1; i<=n_sci_parms; ++i) {
      REAL z = (theta[i] - mean[i])/sdev[i];
      REAL e = minus_log_root_two_pi - log(sdev[i]) - 0.5*pow(z,2);
      sum += den_val(true,e);
    }

    return sum;
  }

  realmat habit_sci_mod::annualize_parms(const realmat& theta)
  {
    habit_sci_mod_parameters mp;
    mp.set_theta(theta);
  
    REAL g = 100.0*12.0*mp.g;
    REAL sig = 100.0*sqrt(12.0)*mp.sig;
    REAL sig_w = 100.0*sqrt(12.0)*mp.sig_w;
    REAL rho = mp.rho;
    REAL phi = pow(mp.phi,12);
    REAL delta = pow(mp.delta,12);
    REAL gamma = mp.gamma;
  
    realmat annual(n_sci_parms,1);
  
    annual[1] = g;
    annual[2] = sig;
    annual[3] = rho;
    annual[4] = sig_w;
    annual[5] = phi;
    annual[6] = delta;
    annual[7] = gamma;
  
    return annual;
  }
  
  realmat habit_sci_mod::annualize_funcs(const realmat& stats)
  {
    REAL  bond_mean  = 100.0*12.0*stats[1];
    REAL  bond_sdev  = 100.0*sqrt(12.0*stats[2]);
    REAL  stock_mean = 100.0*12.0*stats[3];
    REAL  stock_sdev = 100.0*sqrt(12.0*stats[4]);
  
    realmat annual(n_sci_funcs,1);
  
    annual[1] = bond_mean;
    annual[2] = bond_sdev;
    annual[3] = stock_mean;
    annual[4] = stock_sdev;
  
    return annual;
  }

  bool habit_sci_mod::gen_sim(realmat& sim, realmat& stats)
  {
    realmat sim_cgsrpd;
    switch (n_dat_vars) {
      case 1 :
        if (gen_sim_cgsrpd(sim_cgsrpd, stats)) {
          INTEGER N = sim_cgsrpd.ncol();
          sim.resize(1,N); 
          for (INTEGER t=1; t<=N; ++t) {
	    sim(1,t) = sim_cgsrpd(2,t);
          }
          vecwrite("hab.simulation.dat", sim);
          return true;
        }
        return false;
      case 2 :
        if (gen_sim_cgsrpd(sim_cgsrpd, stats)) {
          INTEGER N = sim_cgsrpd.ncol();
          sim.resize(2,N); 
          for (INTEGER t=1; t<=N; ++t) {
	    sim(1,t) = sim_cgsrpd(1,t);
	    sim(2,t) = sim_cgsrpd(2,t);
          }
          vecwrite("hab.simulation.dat", sim);
          return true;
        }
        return false;
      case 3 :
	if (gen_sim_cgsrpd(sim, stats)) {
          vecwrite("hab.simulation.dat", sim);
	  return true;
        }
	return false;
      default :
        return false;
    }
  }
}

