#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 funcs; typedef std::vector::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= 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 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& mod_pfvec, const std::vector& 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& 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