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

Copyright (C) 2018

A. Ronald Gallant
Post Office Box 659
Raleigh 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 "libscl.h"
#include "lprior.h"
#include "slpeqns.h"

using namespace scl;
using namespace std;

int main(int argc, char** argp, char** envp)
{
  const INTEGER L = 1;
  const INTEGER M = 3;
  const INTEGER mrs_pos = 1;
  const INTEGER cf_pos = 3;
  const INTEGER years = 30;
  const INTEGER n_stats = 10;
  const INTEGER n_obs = 5000;
  const REAL lambda = 1.0;

  string filename;
  ifstream fin;

  realmat rho_hat;
  
  if(vecread("dcf.rho_mode.dat", rho_hat)==0) error("Error, vecread failed");

  realmat rho = rho_hat;

  realmat b0, B, S, R, theta;
  frac_rho(M, L, rho, b0, B, S, R, theta);

  realmat E;
  INTEGER ier;
  REAL maxlam = eigen(B,E,ier);

  cout << '\n';
  //cout << "rank = " << rank << '\n';
  cout << "ier = " << ier << '\n';
  cout << "maxlam = " << maxlam << '\n';
  cout << '\n';
  cout << "b0 = " << b0 << '\n';
  cout << "B = " << B << '\n';
  cout << "R = " << R << '\n';
  cout << "S = " << S << '\n';
  cout << "theta = " << theta << '\n';

  REAL beta = theta[1];
  REAL gamma = theta[2];

  realmat I(M,M,0.0);
  for (INTEGER i=1; i<=M; ++i) I(i,i) = 1.0;

  realmat IB = I - B;
  realmat IBinv = inv(IB);

  realmat lim_mean = IBinv*b0;
  realmat vecS(M*M,1);
  for (INTEGER i=1; i<=M*M; ++i) vecS[i] = S[i];
  realmat vecV = inv(kronprd(I,I) - kronprd(B,B))*vecS;
  realmat lim_V(M,M);
  for (INTEGER i=1; i<=M*M; ++i) lim_V[i] = vecV[i];
  realmat lim_R = lim_V;
  if (factor(lim_R)!=0) error("Error, factor failed");

  cout << '\n';
  cout << "lim_mean = " << lim_mean << '\n';
  cout << "lim_V = " << lim_V << '\n';
  cout << "lim_R = " << lim_R << '\n';
  cout << "n_obs " << n_obs << '\n';

  intvec idx(n_stats);
  realmat stats(n_stats,1);
  realmat target(n_stats,1);

  target[ 1] = 0.02;
  target[ 2] = 0.04;
  target[ 3] = 0.01;
  target[ 4] = 0.02;
  target[ 5] = 0.03;
  target[ 6] = 0.02;
  target[ 7] = 0.07;
  target[ 8] = 0.14;
  target[ 9] = 0.07;
  target[10] = 0.04;

  idx[ 1] = 1;
  idx[ 2] = 1;
  idx[ 3] = 1;
  idx[ 4] = 1;
  idx[ 5] = 1;
  idx[ 6] = 1;
  idx[ 7] = 1;
  idx[ 8] = 1;
  idx[ 9] = 1;
  idx[10] = 1;

  denval lp=lprior(M,L,mrs_pos,cf_pos,years,rho,lambda,target,idx,stats);

  cout << '\n';
  cout << boolalpha;
  cout << "lp = (" << lp.positive << ", " << lp.log_den << ")" << '\n';
  cout << '\n';
  cout << "stat        idx    value    target" << '\n';
  cout << "lcg_mean     "  << idx[ 1]
    << fmt('f',10,5,stats[ 1]) << fmt('f',10,5,target[ 1]) << '\n';
  cout << "lcg_sdev     "  << idx[ 2]
    << fmt('f',10,5,stats[ 2]) << fmt('f',10,5,target[ 2]) << '\n';
  cout << "yld_mean_01  "  << idx[ 3]
    << fmt('f',10,5,stats[ 3]) << fmt('f',10,5,target[ 3]) << '\n';
  cout << "yld_sdev_01  "  << idx[ 4]
    << fmt('f',10,5,stats[ 4]) << fmt('f',10,5,target[ 4]) << '\n';
  cout << "yld_mean_yr  "  << idx[ 5]
    << fmt('f',10,5,stats[ 5]) << fmt('f',10,5,target[ 5]) << '\n';
  cout << "yld_sdev_yr  "  << idx[ 6]
    << fmt('f',10,5,stats[ 6]) << fmt('f',10,5,target[ 6]) << '\n';
  cout << "stk_mean_01  "  << idx[ 7]
    << fmt('f',10,5,stats[ 7]) << fmt('f',10,5,target[ 7]) << '\n';
  cout << "stk_sdev_01  "  << idx[ 8]
    << fmt('f',10,5,stats[ 8]) << fmt('f',10,5,target[ 8]) << '\n';
  cout << "stk_mean_yr  "  << idx[ 9]
    << fmt('f',10,5,stats[ 9]) << fmt('f',10,5,target[ 9]) << '\n';
  cout << "stk_sdev_yr  "  << idx[10]
    << fmt('f',10,5,stats[10]) << fmt('f',10,5,target[10]) << '\n';
  cout << '\n';

  REAL implied_gamma = sqrt(lim_V(1,1))/target[2];
  REAL implied_beta = exp(lim_mean[1] + implied_gamma*target[1]);

  cout << "implied gamma " << fmt('f',10,5,implied_gamma) << '\n';
  cout << "implied beta  " << fmt('f',10,5,implied_beta) << '\n';
  cout << '\n';

  realmat y, ecf, pvcf, pv1, cumecf, dcf, yld;

  INT_32BIT seed = 893475;

  realmat z(M,1);
  for (INTEGER i=1; i<=M; ++i) z[i] = unsk(seed);

  realmat ylag = lim_mean + lim_R*z;

  dcfyld(b0,B,S,ylag,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cumecf,dcf,yld);
  REAL price_lag = dcf[years];

  cout << '\n';
  cout << "ylag = " << ylag << '\n';
  cout << "ecf[1] = " << ecf[1] << '\n';
  cout << "ecf[years] = " << ecf[years] << '\n';
  cout << "dcf[1] = " << dcf[1] << '\n';
  cout << "dcf[years] = " << dcf[years] << '\n';
  cout << "yld[1] = " << yld[1] << '\n';
  cout << "yld[years] = " << yld[years] << '\n';

  realmat sim(M,n_obs);
  realmat data(2,n_obs);

  for (INTEGER t=1; t<=n_obs; ++t) {
    for (INTEGER i=1; i<=M; ++i) z[i] = unsk(seed);
    realmat y = b0 + B*ylag + R*z;
    //realmat y = lim_mean + lim_R*z;
    for (INTEGER i=1; i<=M; ++i) sim(i,t) = y[i];
    dcfyld(b0,B,S,y,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cumecf,dcf,yld);
    REAL price = dcf[years];
    data(1,t) = log(price+y[cf_pos]) - log(price_lag);
    data(2,t) = (log(beta) - y[mrs_pos])/gamma;
    ylag = y;
    price_lag = price;
  }
  
  cout << '\n';
  cout << "n_obs " << n_obs << '\n';
  cout << "sim(,1) = " << sim("",1) << '\n';
  cout << "sim(,n_obs) = " << sim("",n_obs) << '\n';
  cout << "data(,1) = " << data("",1) << '\n';
  cout << "data(,n_obs) = " << data("",n_obs) << '\n';

  for (INTEGER t=1; t<=n_obs; ++t) {
    if(!IsFinite(data(1,t))) {
      cout << t << fmt('f',10,5,sim(1,t)) << fmt('f',10,5,sim(2,t)) << '\n';
      cout << t << fmt('f',10,5,data(1,t)) << fmt('f',10,5,data(2,t)) << '\n';
    }
  }

  cout << simple(T(data));


  return 0;
}

