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

Copyright (C) 2018

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.

-------------------------------------------------------------------------------

Function      dcfyld - compute expected cash flows, discounted cash flows, 
                       and yields from varcoef output

Syntax        #include "libscl.h"
Syntax        #include "dfcyld.h"
              void dcfyld(const realmat& b0, const realmat& B, 
                const realmat& S, const realmat& y0, 
                INTEGER mrs_pos, INTEGER cf_pos, INTEGER years, 
		realmat& ecf, realmat& pvcf, realmat& pv1, 
		realmat& cumecf; realmat& dcf, realmat& yld), 

Prototype in  dfcyld.h

Description   Uses as inputs b0, B, S, computed by varcoef, initial 
              condition y0, which is typically y0 = data(t,"") for some t, 
              mrs_pos, that gives the element of y0 containing log mrs, 
              and cf_pos, that gives the element of y0 containing log cash 
              flow growth, and years for which discounted payoffs are 
              desired; dcfyld computes ecf, the expected value of the
	      cash flow for each year, pvcf, the present value of the cash 
              flow for each year, pv1, the present value of a dollar for 
              each year, cumecf, the cumulative expected value of the
	      cash flow, dcf, the cumulative present value of the cash 
	      flow, and yld, the yield curve.

Return value  None

Remark        The timing convention is that, if y0 = data(t,""), pvcf[1] 
              is price at time t of a time t+1 payoff, pvcf[2] is price at 
              time n of a time t+2 payoff, pvcf[3] is price at time t of 
              a time n+t payoff, etc.
              One can handle a VAR with more than one lag by fitting using 
              varcoef with more than one lag and putting the VAR in state
              space form before calling dcfyld.

Reference:    Gallant, A. Ronald, and George Tauchen (2018), "Discounted
              Cash Flows," http://www.aronaldg.org/papers/dcf.pdf

Functions     Library: (none)
called        libscl:  realmat, varcoef

------------------------------------------------------------------------------*/
#include "libscl.h"
#include "dcfyld.h"

using namespace std;
using namespace scl;

void dcfyld(const realmat& b0, const realmat& B, const realmat& S, 
  const realmat& y0, INTEGER mrs_pos, INTEGER cf_pos, INTEGER years, 
  realmat& ecf, realmat& pvcf, realmat& pv1, 
  realmat& cumecf, realmat& dcf, realmat& yld)
{
  const bool debug = false;

  if (debug) {
    cerr << '\n';
    cerr << starbox("/Inputs//");
    cerr << "\t b0 = " << b0 << '\n';
    cerr << "\t B = " << B << '\n';
    cerr << "\t S = " << S << '\n';
    cerr << "\t y0 = " << y0 << '\n';
    cerr << "\t mrs_pos = " << mrs_pos << '\n';
    cerr << "\t cf_pos = " << cf_pos << '\n';
    cerr << '\n';
  }

  INTEGER M = b0.size();

  if (b0.ncol() != 1 && b0.nrow() != 1) error("Error, dcfyld, b0 not a vector");
  if (B.ncol() != M || B.nrow() != M) error("Error, dcfyld, B dimensions bad");
  if (S.nrow() != M || S.ncol() != M) error("Error, dcfyld, S dimensions bad");
  if (y0.size() != M) error("Error, dcfyld, y0 dimensions bad");

  realmat x0(M,1);
  for (INTEGER i=1; i<=M; ++i) x0[i] = y0[i];

  vector<realmat> C(years+1);;
  realmat I(M,M,0.0);
  for (INTEGER i=1; i<=M; ++i) I(i,i) = 1.0;
  C[years] = I;
  for (INTEGER t=1; t<=years; ++t) {
    C[years-t] = I + C[years + 1 - t]*B;
  }

  if (debug) {
    cerr << '\n';
    cerr << starbox("/Coefficient of variance recursion//");
    for (INTEGER t=0; t<=years; ++t) {
      cerr << "\t t = " << t << C[t];
    }
    cerr << '\n';
  }

  vector<realmat> cum_V(years+1);
  cum_V[0] = realmat();
  for (INTEGER t=1; t<=years; ++t) {
    realmat wt_sum_S(M,M,0.0);
    for (INTEGER j=1; j<=t; ++j) {
      wt_sum_S += (C[years-j+1]*S)*T(C[years-j+1]);
    }
    cum_V[t] = wt_sum_S;
  }

  if (debug) {
    cerr << '\n';
    cerr << starbox("/Cumulative variance//");
    for (INTEGER t=0; t<=years; ++t) {
      cerr << "\t t = " << t << cum_V[t];
    }
    cerr << '\n';
  }

  vector<realmat> b(years+1);
  vector<realmat> cum_b(years+1);
  b[0] = realmat();
  b[1] = b0;
  cum_b[0] = realmat();
  cum_b[1] = b0;
  for (INTEGER t=2; t<=years; ++t) {
    b[t] = b0 + B*b[t-1];
    cum_b[t] = cum_b[t-1] + b[t];
  }

  if (debug) {
    cerr << '\n';
    cerr << starbox("/b recursion//");
    for (INTEGER t=0; t<=years; ++t) {
      cerr << "\t t = " << t << b[t];
    }
    cerr << '\n';
  }

  if (debug) {
    cerr << '\n';
    cerr << starbox("/cum_b//");
    for (INTEGER t=0; t<=years; ++t) {
      cerr << "\t t = " << t << cum_b[t];
    }
    cerr << '\n';
  }

  vector<realmat> y(years+1);;
  vector<realmat> cum_y(years+1);;
  y[0] = x0;
  cum_y[0] = realmat(M,1,0.0);
  for (INTEGER t=1; t<=years; ++t) {
    y[t] = B*y[t-1];
    cum_y[t] = cum_y[t-1] + y[t];
  }

  if (debug) {
    cerr << '\n';
    cerr << starbox("/y recursion//");
    for (INTEGER t=0; t<=years; ++t) {
      cerr << "\t t = " << t << y[t];
    }
    cerr << '\n';
  }

  if (debug) {
    cerr << '\n';
    cerr << starbox("/cum_y//");
    for (INTEGER t=0; t<=years; ++t) {
      cerr << "\t t = " << t << cum_y[t];
    }
    cerr << '\n';
  }

  realmat M_mean(years,1);
  //realmat M_var(years,1);
  //realmat M_sdev(years,1);

  for (INTEGER t=1; t<=years; ++t) {
    realmat M = cum_b[t] + cum_y[t];
    realmat V = cum_V[t];
    REAL mean = M[mrs_pos];
    REAL var = V(mrs_pos,mrs_pos);
    M_mean[t] = exp(mean + 0.5*var);
    //M_var[t] = exp(2.0*mean + 2.0*var) - exp(2.0*mean + var);
    //M_sdev[t] = sqrt(var);
  }

  pv1.resize(years,1);
  yld.resize(years,1);
  for (INTEGER t=1; t<=years; ++t) {
    pv1[t] = M_mean[t];
    yld[t] = -log(M_mean[t])/REAL(t);
  }

  realmat dcf_mean(years,1);
  //realmat dcf_var(years,1);
  //realmat dcf_sdev(years,1);

  for (INTEGER t=1; t<=years; ++t) {
    realmat M = cum_b[t] + cum_y[t];
    realmat V = cum_V[t];
    REAL mean = M[mrs_pos] + M[cf_pos];
    REAL var = V(mrs_pos,mrs_pos) + 2.0*V(mrs_pos,cf_pos) + V(cf_pos,cf_pos);
    dcf_mean[t] = exp(mean + 0.5*var);
    //dcf_var[t] = exp(2.0*mean + 2.0*var) - exp(2.0*mean + var);
    //dcf_sdev[t] = sqrt(var);
  }

  pvcf.resize(years,1);
  dcf.resize(years,1);
  REAL sum = 0.0;
  for (INTEGER t=1; t<=years; ++t) {
    pvcf[t] = dcf_mean[t];
    sum += dcf_mean[t];
    dcf[t] = sum;
  }

  realmat ecf_mean(years,1);
  //realmat ecf_var(years,1);
  //realmat ecf_sdev(years,1);

  for (INTEGER t=1; t<=years; ++t) {
    realmat M = cum_b[t] + cum_y[t];
    realmat V = cum_V[t];
    REAL mean = M[cf_pos];
    REAL var = V(cf_pos,cf_pos);
    ecf_mean[t] = exp(mean + 0.5*var);
    //ecf_var[t] = exp(2.0*mean + 2.0*var) - exp(2.0*mean + var);
    //ecf_sdev[t] = sqrt(var);
  }

  cumecf.resize(years,1);
  ecf.resize(years,1);
  sum = 0.0;
  for (INTEGER t=1; t<=years; ++t) {
    ecf[t] = ecf_mean[t];
    sum += ecf_mean[t];
    cumecf[t] = sum;
  }

  return;
}
