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

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> 
  lrr_sci_mod_variables::get_lrr_sci_mod_variables(realmat& mv)
  {
    vector<string> names;
    names.push_back("log_consumption_growth"); mv = log_consumption_growth;
    names.push_back("consumption_growth"); mv = cbind(mv,consumption_growth);
    names.push_back("consumption"); mv = cbind(mv,consumption);
    names.push_back("long_run_risk"); mv = cbind(mv,long_run_risk);
    names.push_back("squared_volatility"); mv = cbind(mv,squared_volatility);
    names.push_back("volatility"); mv = cbind(mv,volatility);
    names.push_back("log_dividend_growth"); mv = cbind(mv,log_dividend_growth);
    names.push_back("dividend_growth"); mv = cbind(mv,dividend_growth);
    names.push_back("dividend"); mv = cbind(mv,dividend);
    names.push_back("zeta[,1]"); 
      names.push_back("zeta[,2]"); 
      names.push_back("zeta[,3]"); 
      mv = cbind(mv,zeta);
    names.push_back("Y[,1]"); 
      names.push_back("Y[,2]"); 
      names.push_back("Y[,3]"); 
      mv = cbind(mv,Y);
    names.push_back("log_mrs"); mv = cbind(mv,log_mrs);
    names.push_back("mrs"); mv = cbind(mv,mrs);
    names.push_back("log_pc_ratio"); mv = cbind(mv,log_pc_ratio);
    names.push_back("geometric_consumption_return"); 
      mv = cbind(mv,geometric_consumption_return);
    names.push_back("geometric_risk_free_rate"); 
      mv = cbind(mv,geometric_risk_free_rate);
    names.push_back("log_pd_ratio"); mv = cbind(mv,log_pd_ratio);
    names.push_back("geometric_stock_return"); 
      mv = cbind(mv,geometric_stock_return);
    return names;
  }

  lrr_sci_mod::lrr_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 = lrr_sci_mod_variables(n0,n,n1);

    #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();

      realmat theta;
      mp.get_theta(theta);
      detail << starbox("/theta//") << '\n';
      detail << mp;
      detail << starbox("/implied annual values//");
      realmat annual = annualize_parms(theta);
      detail << cbind(theta,annual);
      detail << '\n';
      detail.flush();

      detail << starbox("/begin gen_sim_cgsr output//") << '\n';;

      realmat sim, stats;
      gen_sim_cgsr(sim,stats);
      vecwrite("lrr.simulation.dat",sim);

      detail << starbox("/end gen_sim_cgsr output//") << '\n';;
      detail.flush();

      realmat variables;
      vector<string> names = mv.get_lrr_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',17,8,sqrt(var[j])) << '\n';
         }
         else {
          detail << '\t' << names[j-1] << ' ' 
            << fmt('f',45-names[j-1].size(),5,mean[j]) 
            << fmt('f',17,8,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
  }
  
  ostream& operator<<(ostream& os, const lrr_sci_mod_parameters& mp)
  {
    os << "\t 1 delta  = " << scl::fmt('f',12,8,mp.delta) << '\n';
    os << "\t 2 gamma  = " << scl::fmt('f',12,8,mp.gamma) << '\n';
    os << "\t 3 psi    = " << scl::fmt('f',12,8,mp.psi) << '\n';
    os << "\t 4 mu_c   = " << scl::fmt('f',12,8,mp.mu_c) << '\n';
    os << "\t 5 rho    = " << scl::fmt('f',12,8,mp.rho) << '\n';
    os << "\t 6 phi_e  = " << scl::fmt('f',12,8,mp.phi_e) << '\n';
    os << "\t 7 varbar = " << scl::fmt('f',12,8,mp.varbar) << '\n';
    os << "\t 8 nu     = " << scl::fmt('f',12,8,mp.nu) << '\n';
    os << "\t 9 sig_w  = " << scl::fmt('f',12,8,mp.sig_w) << '\n';
    os << "\t10 mu_d   = " << scl::fmt('f',12,8,mp.mu_d) << '\n';
    os << "\t11 phi_d  = " << scl::fmt('f',12,8,mp.phi_d) << '\n';
    os << "\t12 pi_d   = " << scl::fmt('f',12,8,mp.pi_d) << '\n';
    os << "\t13 phi_u  = " << scl::fmt('f',12,8,mp.phi_u) << '\n';
    os << '\n';
    os << "\t   rpsi   = " << scl::fmt('f',12,8,mp.rpsi) << '\n';
    os << "\t   theta  = " << scl::fmt('f',12,8,mp.theta) << '\n';
    return os;
  }

  bool lrr_sci_mod::gen_sim_cgsr(realmat& sim)
  {
  
    INT_32BIT seed = 47;
  
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << starbox("/model variables and parameters//") << '\n';
      debug << "\t mv.n0 = " << mv.n0 << '\n';
      debug << "\t mv.n  = " << mv.n  << '\n';
      debug << "\t mv.n1 = " << mv.n1 << '\n';
      debug << "\t mv.nt = " << mv.nt << '\n';
      debug << '\n';
      debug << mp << '\n';
      debug << starbox("/means and std. dev. of model variables//") << '\n';
      debug.flush();
    }
    #endif
  
    REAL C = REAL_MIN;
    REAL D = REAL_MIN;
    for (INTEGER t=1; t<mv.nt; ++t) {
  
      REAL var_t = mv.squared_volatility[t];
      REAL sig_t = mv.volatility[t];
      REAL x_t   = mv.long_run_risk[t];
  
      mv.Y(t,1) = 1.0;
      mv.Y(t,2) = x_t;
      mv.Y(t,3) = var_t;
  
      REAL eta = unsk(seed);
  
      REAL del_c = mp.mu_c + x_t + sig_t*eta;
      REAL cg = exp(del_c);
      C *= cg;
  
      mv.log_consumption_growth[t+1] = del_c;
      mv.consumption_growth[t+1] = cg; 
      mv.consumption[t+1] = C;
  
      REAL e = unsk(seed);
      mv.long_run_risk[t+1] = mp.rho*x_t + mp.phi_e*sig_t*e;
  
      REAL w = unsk(seed);
  
      REAL var = mp.varbar + mp.nu*(var_t - mp.varbar) + mp.sig_w*w;
      REAL sig = sqrt(var > 0.0 ? var : 0.0);
  
      mv.zeta(t+1,1) = sig_t*eta;
      mv.zeta(t+1,2) = sig_t*e;
      mv.zeta(t+1,3) = mp.sig_w*w;
  
      mv.squared_volatility[t+1] = var;
      mv.volatility[t+1] = sig;
  
      REAL u = unsk(seed);
  
      REAL del_d = mp.mu_d + mp.phi_d*x_t + mp.pi_d*sig_t*eta + mp.phi_u*sig_t*u;
      REAL dg = exp(del_d);
      D *= dg;
  
      mv.log_dividend_growth[t+1] = del_d;
      mv.dividend_growth[t+1] = dg;
      mv.dividend[t+1] = D;
    }
    if (!IsFinite(C)) error("Error, sci_mod::make_state, C not finite");
    if (!IsFinite(D)) error("Error, sci_mod::make_state, D not finite");
  
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      REAL lcg_sum = 0.0;
      REAL cg_sum = 0.0;
      REAL x_sum = 0.0;
      REAL sv_sum = 0.0;
      REAL ldg_sum = 0.0;
      REAL dg_sum = 0.0;
      for (INTEGER t=mv.n0+1; t<=mv.n0+mv.n; ++t) {
        lcg_sum += mv.log_consumption_growth[t];
        cg_sum += mv.consumption_growth[t];
        x_sum += mv.long_run_risk[t];
        sv_sum += mv.squared_volatility[t];
        ldg_sum += mv.log_dividend_growth[t];
        dg_sum += mv.dividend_growth[t];
      }
      REAL lcg_mean = lcg_sum/REAL(mv.n);
      REAL cg_mean = cg_sum/REAL(mv.n);
      REAL x_mean = x_sum/REAL(mv.n);
      REAL sv_mean = sv_sum/REAL(mv.n);
      REAL ldg_mean = ldg_sum/REAL(mv.n);
      REAL dg_mean = dg_sum/REAL(mv.n);
  
      lcg_sum = 0.0;
      cg_sum = 0.0;
      x_sum = 0.0;
      sv_sum = 0.0;
      ldg_sum = 0.0;
      dg_sum = 0.0;
      for (INTEGER t=mv.n0+1; t<=mv.n0+mv.n; ++t) {
        lcg_sum += pow(mv.log_consumption_growth[t]-lcg_mean,2);
        cg_sum += pow(mv.consumption_growth[t]-cg_mean,2);
        x_sum += pow(mv.long_run_risk[t]-x_mean,2);
        sv_sum += pow(mv.squared_volatility[t]-sv_mean,2);
        ldg_sum += pow(mv.log_dividend_growth[t]-ldg_mean,2);
        dg_sum += pow(mv.dividend_growth[t]-dg_mean,2);
      }
      REAL lcg_sdev = sqrt(lcg_sum/REAL(mv.n-1));
      REAL cg_sdev = sqrt(cg_sum/REAL(mv.n-1));
      REAL x_sdev = sqrt(x_sum/REAL(mv.n-1));
      REAL sv_sdev = sqrt(sv_sum/REAL(mv.n-1));
      REAL ldg_sdev = sqrt(ldg_sum/REAL(mv.n-1));
      REAL dg_sdev = sqrt(dg_sum/REAL(mv.n-1));
  
      debug << "\t lcg_mean = " << fmt('f',12,8,lcg_mean) 
           << ",  lcg_sdev = " << fmt('f',12,8,lcg_sdev) << '\n';
      debug << "\t cg_mean =  " << fmt('f',12,8,cg_mean) 
           << ",  cg_sdev =  " << fmt('f',12,8,cg_sdev) << '\n';
      debug << "\t x_mean =   " << fmt('f',12,8,x_mean) 
           << ",  x_sdev =   " << fmt('f',12,8,x_sdev) << '\n';
      debug << "\t sv_mean =  " << fmt('f',12,8,sv_mean) 
           << ",  sv_sdev =  " << fmt('f',12,8,sv_sdev) << '\n';
      debug << "\t ldg_mean = " << fmt('f',12,8,ldg_mean) 
           << ",  ldg_sdev = " << fmt('f',12,8,ldg_sdev) << '\n';
      debug << "\t dg_mean =  " << fmt('f',12,8,dg_mean) 
           << ",  dg_sdev  = " << fmt('f',12,8,dg_sdev) << '\n';
      debug.flush();
    }
    #endif
  
    // convergence tolerance
    const REAL tol = 1.0e-6;
  
    // temporaries
    REAL z, A2_1, A2_2, A0_1, A0_2, A0_3, A0_4;

    // consumption contraction mapping
    INTEGER count = 0;
    REAL diff = REAL_MAX;
    REAL z_bar = 7.0;
    REAL k1, k0;
    REAL A0, A1, A2;
    while (diff > tol) {
      ++count;
      if (count > 5000) return false;
      z = z_bar;
      k1 = exp(z)/(1.0 + exp(z));
      k0 = log(1.0 + exp(z)) - k1*z;
      A2_1 = -(mp.gamma - 1.0)*(1.0 - mp.rpsi)/(2.0*(1.0 - k1*mp.nu));
      A2_2 = 1.0 + pow(k1*mp.phi_e/(1.0 - k1*mp.rho),2);
      A2 = A2_1*A2_2;
      A1 = (1.0 - mp.rpsi)/(1.0 - k1*mp.rho);
      A0_1 = 1.0/(1.0 - k1);
      A0_2 = log(mp.delta) + k0 + (1.0 - mp.rpsi)*mp.mu_c;
      A0_3 = k1*A2*(1.0 - mp.nu)*mp.varbar;
      A0_4 = (mp.theta/2.0)*pow(k1*A2*mp.sig_w,2);
      A0 = A0_1*(A0_2 + A0_3 + A0_4);
      z_bar = A0 + mp.varbar*A2;
      diff = fabs(z_bar - z);
    }
    z_bar = z;
  
  
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << "\t count =    " << fmt('d',12,count) 
           << ",  z_bar =    " << fmt('f',12,8,z_bar) << '\n';
    }
    #endif 
  
    for (INTEGER t=1; t<=mv.nt; ++t) {
      mv.log_pc_ratio[t] 
        = A0 + A1*mv.long_run_risk[t] + A2*mv.squared_volatility[t];
    }
  
    // parameters for mrs
    REAL lambda_eta, lambda_e, lambda_w;
    REAL lambda_w_1, lambda_w_2;
    REAL Gamma0, Gamma1, Gamma2;
    REAL Gamma0_1, Gamma0_2;
    Gamma0_1 = log(mp.delta) - mp.rpsi*mp.mu_c;
    Gamma0_2 = -0.5*mp.theta*(mp.theta - 1.0)*pow(k1*A2*mp.sig_w,2);
    Gamma0 = Gamma0_1 + Gamma0_2;
    Gamma1 = -mp.rpsi;
    Gamma2 = (mp.theta - 1.0)*(k1*mp.nu - 1.0)*A2;
    lambda_eta = mp.gamma;
    lambda_e = (mp.gamma - mp.rpsi)*(k1*mp.phi_e)/(1.0 - k1*mp.rho);
    lambda_w_1 = -(mp.gamma - 1.0)*(mp.gamma - mp.rpsi)*0.5*k1/(1.0 - k1*mp.nu);
    lambda_w_2 = 1.0 + pow(k1*mp.phi_e/(1.0 - k1*mp.rho),2);
    lambda_w = lambda_w_1*lambda_w_2;
  
    // mrs, geometric consumption return
    for (INTEGER t=1; t<mv.nt; ++t) {
      mv.log_mrs[t+1] 
        = Gamma0 + Gamma1*mv.Y(t,2) + Gamma2*mv.Y(t,3)
          - lambda_eta*mv.zeta(t+1,1) - lambda_e*mv.zeta(t+1,2)
          - lambda_w*mv.zeta(t+1,3);
      mv.mrs[t+1] = exp(mv.log_mrs[t+1]);
      mv.geometric_consumption_return[t+1] 
        = k0 + k1*mv.log_pc_ratio[t+1] + mv.log_consumption_growth[t+1]
          - mv.log_pc_ratio[t];
    }
  
    // geometric risk free rate
    for (INTEGER t=1; t<=mv.nt; ++t) {
      mv.geometric_risk_free_rate[t]
        = -Gamma0 - Gamma1*mv.Y(t,2) - Gamma2*mv.Y(t,3)
          - 0.5*pow(lambda_eta,2)*mv.squared_volatility[t]
          - 0.5*pow(lambda_e,2)*mv.squared_volatility[t]
          - 0.5*pow(lambda_w*mp.sig_w,2);
    }
  
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      REAL mrs_sum = 0.0;
      REAL lpc_sum = 0.0;
      REAL lcr_sum = 0.0;
      REAL lrf_sum = 0.0;
      for (INTEGER t=mv.n0+1; t<=mv.n0+mv.n; ++t) {
        mrs_sum += mv.mrs[t];
        lpc_sum += mv.log_pc_ratio[t];
        lcr_sum += mv.geometric_consumption_return[t];
        lrf_sum += mv.geometric_risk_free_rate[t];
      }
      REAL mrs_mean = mrs_sum/REAL(mv.n);
      REAL lpc_mean = lpc_sum/REAL(mv.n);
      REAL lcr_mean = lcr_sum/REAL(mv.n);
      REAL lrf_mean = lrf_sum/REAL(mv.n);
    
      mrs_sum = 0.0;
      lpc_sum = 0.0;
      lcr_sum = 0.0;
      lrf_sum = 0.0;
      for (INTEGER t=mv.n0+1; t<=mv.n0+mv.n; ++t) {
        mrs_sum += pow(mv.mrs[t]-mrs_mean,2);
        lpc_sum += pow(mv.log_pc_ratio[t]-lpc_mean,2);
        lcr_sum += pow(mv.geometric_consumption_return[t]-lcr_mean,2);
        lrf_sum += pow(mv.geometric_risk_free_rate[t]-lrf_mean,2);
      }
      REAL mrs_sdev = sqrt(mrs_sum/REAL(mv.n-1));
      REAL lpc_sdev = sqrt(lpc_sum/REAL(mv.n-1));
      REAL lcr_sdev = sqrt(lcr_sum/REAL(mv.n-1));
      REAL lrf_sdev = sqrt(lrf_sum/REAL(mv.n-1));
    
      debug << "\t mrs_mean = " << fmt('f',12,8,mrs_mean) 
           << ",  mrs_sdev = " << fmt('f',12,8,mrs_sdev) << '\n';
      debug << "\t lpc_mean = " << fmt('f',12,8,lpc_mean) 
           << ",  lpc_sdev = " << fmt('f',12,8,lpc_sdev) << '\n';
      debug << "\t lcr_mean = " << fmt('f',12,8,lcr_mean) 
           << ",  lcr_sdev = " << fmt('f',12,8,lcr_sdev) << '\n';
      debug << "\t lrf_mean = " << fmt('f',12,8,lrf_mean) 
           << ",  lrf_sdev = " << fmt('f',12,8,lrf_sdev) << '\n';
      debug.flush();
    }
    #endif
  
    // dividend contraction mapping
    count = 0;
    diff = REAL_MAX;
    z_bar = 5.0;
    REAL k0d, k1d;
    REAL A0d, A1d, A2d;
    while (diff > tol) {
      ++count;
      if (count > 5000) return false;
      z = z_bar;
      k1d = exp(z)/(1.0 + exp(z));
      k0d = log(1.0 + exp(z)) - k1d*z;
      A1d = (mp.phi_d - mp.rpsi)/(1.0 - k1d*mp.rho);
      A2_1 = 1.0/(1.0 - k1d*mp.nu);
      A2_2 = Gamma2 + 0.5*pow(mp.pi_d-lambda_eta,2) 
             + 0.5*pow(k1d*A1d*mp.phi_e-lambda_e,2);
      A2d = A2_1*A2_2;
      A0_1 = 1.0/(1.0 - k1d);
      A0_2 = Gamma0 + k0d + mp.mu_d + k1d*A2d*(1.0 - mp.nu)*mp.varbar;
      A0_3 = 0.5*pow(k1d*A2d - lambda_w,2)*pow(mp.sig_w,2);
      A0d = A0_1*(A0_2 + A0_3);
      z_bar = A0d + mp.varbar*A2d;
      diff = fabs(z_bar - z);
    }
    z_bar = z;
  
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      debug << "\t count =    " << fmt('d',12,count) 
           << ",  z_bar =    " << fmt('f',12,8,z_bar) << '\n';
    }
    #endif
  
    for (INTEGER t=1; t<=mv.nt; ++t) {
      mv.log_pd_ratio[t] 
        = A0d + A1d*mv.long_run_risk[t] + A2d*mv.squared_volatility[t];
    }
  
    // geometric dividend return
    for (INTEGER t=1; t<mv.nt; ++t) {
      mv.geometric_stock_return[t+1] 
        = k0d + k1d*mv.log_pd_ratio[t+1] + mv.log_dividend_growth[t+1]
          - mv.log_pd_ratio[t];
    }
  
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      REAL lpd_sum = 0.0;
      REAL ldr_sum = 0.0;
      for (INTEGER t=mv.n0+1; t<=mv.n0+mv.n; ++t) {
        lpd_sum += mv.log_pd_ratio[t];
        ldr_sum += mv.geometric_stock_return[t];
      }
      REAL lpd_mean = lpd_sum/REAL(mv.n);
      REAL ldr_mean = ldr_sum/REAL(mv.n);
    
      lpd_sum = 0.0;
      ldr_sum = 0.0;
      for (INTEGER t=mv.n0+1; t<=mv.n0+mv.n; ++t) {
        lpd_sum += pow(mv.log_pd_ratio[t]-lpd_mean,2);
        ldr_sum += pow(mv.geometric_stock_return[t]-ldr_mean,2);
      }
      REAL lpd_sdev = sqrt(lpd_sum/REAL(mv.n-1));
      REAL ldr_sdev = sqrt(ldr_sum/REAL(mv.n-1));
    
      debug << "\t lpd_mean = " << fmt('f',12,8,lpd_mean) 
           << ",  lpd_sdev = " << fmt('f',12,8,lpd_sdev) << '\n';
      debug << "\t ldr_mean = " << fmt('f',12,8,ldr_mean) 
           << ",  ldr_sdev = " << fmt('f',12,8,ldr_sdev) << '\n';
      debug.flush();
    }
    #endif
  
    INTEGER N = mv.n/12;
  
    sim.resize(2,N);
  
    REAL r0, C0, C12;
    INTEGER tt = 0;
    for (INTEGER t=1; t <= N; ++t) {
      r0 = C0 = C12 = 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];
      }
      sim(1,t) = log(C0) - log(C12);
      sim(2,t) = r0;
    }
    
    #if defined DEBUG_GMUSR
    if (debug_switch) {
      realmat ladg(N,1);
      realmat larf(N,1);
      REAL D0, D12;
      tt = 0;
      for (INTEGER t=1; t <= N; ++t) {
        r0 = D0 = D12 = 0.0;
        for (INTEGER j=1; j <= 12; j++) {
          ++tt;
          r0  += mv.geometric_risk_free_rate[mv.n0 + tt];
          D0  += mv.dividend[mv.n0 + tt];
          D12 += mv.dividend[mv.n0 + tt - 12];
        }
        ladg[t] = log(D0) - log(D12);
        larf[t] = r0;
      }
    
      REAL lacg_sum = 0.0;
      REAL ladg_sum = 0.0;
      REAL ladr_sum = 0.0;
      REAL larf_sum = 0.0;
      for (INTEGER t=1; t<=N; ++t) {
        lacg_sum += sim(1,t);
        ladg_sum += ladg[t];
        ladr_sum += sim(2,t);
        larf_sum += larf[t];
      }
      REAL lacg_mean = lacg_sum/REAL(N);
      REAL ladg_mean = ladg_sum/REAL(N);
      REAL ladr_mean = ladr_sum/REAL(N);
      REAL larf_mean = larf_sum/REAL(N);
    
      lacg_sum = 0.0;
      ladg_sum = 0.0;
      ladr_sum = 0.0;
      larf_sum = 0.0;
      for (INTEGER t=1; t<=N; ++t) {
        lacg_sum += pow(sim(1,t)-lacg_mean,2);
        ladg_sum += pow(ladg[t]-ladg_mean,2);
        ladr_sum += pow(sim(2,t)-ladr_mean,2);
        larf_sum += pow(larf[t]-larf_mean,2);
      }
      REAL lacg_sdev = sqrt(lacg_sum/REAL(N-1));
      REAL ladg_sdev = sqrt(ladg_sum/REAL(N-1));
      REAL ladr_sdev = sqrt(ladr_sum/REAL(N-1));
      REAL larf_sdev = sqrt(larf_sum/REAL(N-1));
    
      debug << "\t lacg_mean =" << fmt('f',12,8,lacg_mean) 
           << ",  lacg_sdev =" << fmt('f',12,8,lacg_sdev) << '\n';
      debug << "\t ladg_mean =" << fmt('f',12,8,ladg_mean) 
           << ",  ladg_sdev =" << fmt('f',12,8,ladg_sdev) << '\n';
      debug << "\t ladr_mean =" << fmt('f',12,8,ladr_mean) 
           << ",  ladr_sdev =" << fmt('f',12,8,ladr_sdev) << '\n';
      debug << "\t larf_mean =" << fmt('f',12,8,larf_mean) 
           << ",  larf_sdev =" << fmt('f',12,8,larf_sdev) << '\n';
      debug.flush();
     
    }
    #endif

    return true;
  }

  bool lrr_sci_mod::gen_sim_cgsr(realmat& sim, realmat& stats)
  {
    bool success = this->gen_sim_cgsr(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 lrr_sci_mod::support(const realmat& theta) 
  {

    realmat theta_lo(n_sci_parms,1);
    
    theta_lo[1]  =  0.9;      // delta
    theta_lo[2]  =  0.25;     // gamma
    theta_lo[3]  =  1.0001;   // psi
    theta_lo[4]  = -0.02;     // mu_c
    theta_lo[5]  =  0.0;      // rho
    theta_lo[6]  =  0.0;      // phi_e
    theta_lo[7]  =  0.0;      // varbar
    theta_lo[8]  =  0.0;      // nu
    theta_lo[9]  =  0.0;      // sig_w
    theta_lo[10] = -0.02;     // mu_d
    theta_lo[11] =  0.0;      // phi_d
    theta_lo[12] =  0.0;      // pi_d
    theta_lo[13] =  0.0;      // pi_u
  
    realmat theta_hi(n_sci_parms,1);
  
    theta_hi[1] =  0.99999;   // delta
    theta_hi[2] =  100.0;     // gamma
    theta_hi[3] =  5.0;       // psi
    theta_hi[4] =  0.02;      // mu_c
    theta_hi[5] =  1.0;       // rho
    theta_hi[6] =  1.0;       // phi_e
    theta_hi[7] =  1.0;       // varbar
    theta_hi[8] =  0.9995;    // nu
    theta_hi[9] =  1.0;       // sig_w
    theta_hi[10] = 0.02;      // mu_d
    theta_hi[11] = 100.0;     // phi_d
    theta_hi[12] = 100.0;     // pi_d
    theta_hi[13] = 100.0;     // pi_u
  
    return ( (theta_lo < theta) && (theta < theta_hi) );
  }

  den_val lrr_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.9992;        // delta
    mean[ 2] =  10.0;           // gamma   
    mean[ 3] =   1.5;           // psi     
    mean[ 4] =   0.0015;        // mu_c    
    mean[ 5] =   0.983;         // rho     
    mean[ 6] =   0.032;         // phi_e   
    mean[ 7] =   0.00004096;    // varbar  
    mean[ 8] =   0.988;         // nu      
    mean[ 9] =   0.0000017;     // sig_w     // Kiku market   
    mean[10] =   0.0012;        // mu_d    
    mean[11] =   2.8;           // phi_d   
    mean[12] =   4.125;         // pi_d    
    mean[13] =   6.2637;        // phi_u   

    /*
    mean[ 9] =   0.0000040;     // sig_w     // Kiku supplied 
    mean[10] =   0.0015;        // mu_d    
    mean[11] =   3.0;           // phi_d   
    mean[12] =   1.8;           // pi_d    
    mean[13] =   4.1243;        // phi_u   
    */

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

    sdev[5] = fabs(mean[5])*0.01/1.96;     // Don't forget to report
    sdev[8] = fabs(mean[8])*0.01/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);
    }

    const REAL b = 0.995;
    const REAL s = 1000.0;
    const REAL logc = log(1.0/s + b);

    REAL nu = theta[8];                      // Don't forget to report
    if (nu < b) {                            // this in the paper.
      sum += den_val(true,-logc);
    }
    else {
      sum += den_val(true,s*(b-nu) - logc);
    }
  
    REAL rho = theta[5];                     // Don't forget to report
    if (rho < b) {                           // this in the paper.
      sum += den_val(true,-logc);
    }
    else {
      sum += den_val(true,s*(b-rho) - logc);
    }
    return sum;
  }

  realmat lrr_sci_mod::annualize_parms(const realmat& theta)
  {
    lrr_sci_mod_parameters mp;
    mp.set_theta(theta);
  
    REAL delta = pow(mp.delta,12);
    REAL gamma = mp.gamma;
    REAL psi = mp.psi;
    REAL mu_c = 100.0*12.0*mp.mu_c;
    REAL rho = pow(mp.rho,12);
    REAL phi_e = mp.phi_e;
    REAL varbar = 12.0*mp.varbar;
    REAL nu = pow(mp.nu,12);
    REAL sig_w = 100.0*sqrt(12.0)*mp.sig_w;
    REAL mu_d = 100.0*12.0*mp.mu_c;
    REAL phi_d = mp.phi_d;
    REAL pi_d = mp.pi_d;
    REAL phi_u = mp.phi_u;

    realmat annual(n_sci_parms,1,0.0);
  
    annual[1]  = delta;
    annual[2]  = gamma;
    annual[3]  = psi;
    annual[4]  = mu_c;
    annual[5]  = rho;
    annual[6]  = phi_e;
    annual[7]  = varbar;
    annual[8]  = nu;
    annual[9]  = sig_w;
    annual[10] = mu_d;
    annual[11] = phi_d;
    annual[12] = pi_d;
    annual[13] = phi_u;
  
    return annual;
  }
  
  realmat lrr_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 lrr_sci_mod::gen_sim(realmat& sim, realmat& stats)
  {
    realmat sim_cgsr;
    switch (n_dat_vars) {
      case 1 :
	if (gen_sim_cgsr(sim_cgsr, stats)) {
	  INTEGER N = sim_cgsr.ncol();
	  sim.resize(1,N); 
	  for (INTEGER t=1; t<=N; ++t) sim[t] = sim_cgsr(2,t);
          return true;
	}
        return false;
      case 2 :
        return  gen_sim_cgsr(sim, stats);
      default :
        return false;
    }
  }

}