#include "libscl.h"
#include "spline_interpolator.h"
using namespace scl;
using namespace std;

extern bool tree_policy(REAL alpha, REAL sigma, REAL beta, REAL gamma,
  spline_interpolator& P, spline_interpolator& Q,
  spline_interpolator& p, spline_interpolator& q);

int main(int argc, char** argp, char** envp)
{
  /*
  REAL alpha = 0.9;
  REAL sigma = 0.1;
  REAL beta  = 0.95;
  REAL gamma = 2.0;

  REAL alpha = 0.95;
  REAL sigma = 0.03;
  REAL beta  = 0.95;
  REAL r = -0.02015;
  REAL gamma = sqrt(-2.0*(r + log(beta)))/sigma;  // Gives risk free rate of r
  */

  REAL alpha = 0.95;
  REAL sigma = 0.02;
  REAL beta  = 0.95;
  REAL gamma = 12.5;
  
  spline_interpolator P;
  spline_interpolator Q;
  spline_interpolator p;
  spline_interpolator q;

  bool converge = tree_policy(alpha, sigma, beta, gamma, P, Q, p, q);

  cout << "converge = " << boolalpha << converge << '\n';

  cout << "linear appoximation to p(y)" << '\n';
  cout << "p = " << p(-0.5) << " + " << p(.5) - p(-.5) << "*(y + 0.5)" << '\n';

  INTEGER N = P.size();
  realmat grid(N,1);
  realmat values(N,1);
  realmat log_values(N,1);
  realmat table(N,4);
  vector<string> header(4);

  header[0] = "y";
  header[1] = "P";
  header[2] = "C";
  header[3] = "p";

  for (INTEGER i=0; i<N; ++i) {
    table(i+1,1) = grid[i+1] = P[i].origin();
    table(i+1,2) = values[i+1] = P(grid[i+1]);
    table(i+1,3) = exp(grid[i+1]);
    table(i+1,4) = log_values[i+1] = p(grid[i+1]);
  }

  realmat errors(2*N,4,0.0);
  for (INTEGER i=1; i<N; ++i) {
    REAL y0 = grid[i];
    REAL y1 = grid[i+1];
    REAL Y1 = exp(y1);
    REAL CG = exp(y1 - y0);
    errors(i,1) = P(y0);
    errors(i,2) = CG;
    errors(i,3) = P(y0) - beta*pow(CG,-gamma)*(Y1 + P(y1));
    errors(i,4) = errors(i,3)/errors(i,1);
  }
  for (INTEGER i=1; i<N; ++i) {
    REAL y1 = grid[i];
    REAL y0 = grid[i+1];
    REAL Y1 = exp(y1);
    REAL CG = exp(y1 - y0);
    errors(N+i,1) = P(y0);
    errors(N+i,2) = CG;
    errors(N+i,3) = P(y0) - beta*pow(CG,-gamma)*(Y1 + P(y1));
    errors(N+i,4) = errors(N+i,3)/errors(N+i,1);
  }

  cout << "Pricing errors " << errors << simple(errors) << '\n';

  /*
  for (INTEGER i=1; i<=N; ++i) {
    cout << "grid[" << fmt('d',2,i) <<"] =" << fmt('f',8,3,grid[i]) << ";  "
         << "values[" << fmt('d',2,i) <<"] =" 
	 << fmt('f',8,3,log_values[i]) <<";\n";
  }
  */

  writetable("P.txt",table,header,15,10);

  header[0] = "P";
  header[1] = "y";
  header[2] = "C";
  header[3] = "y_alt";

  for (INTEGER i=0; i<N; ++i) {
    table(i+1,1) = grid[i+1] = Q[i].origin();
    table(i+1,2) = values[i+1] = Q(grid[i+1]);
    table(i+1,3) = exp(values[i+1]);
    table(i+1,4) = log_values[i+1] = q(log(grid[i+1]));
  }

  writetable("Q.txt",table,header,15,10);

  INTEGER n = 100;

  table.resize(n,4);

  header[0] = "y";
  header[1] = "P";
  header[2] = "C";
  header[3] = "p";

  REAL lo = -4.0*sigma/sqrt(1.0 - alpha*alpha);
  REAL hi = -lo;


  for (INTEGER i=1; i<=n; ++i) {
    table(i,1) = lo + (REAL(i)/REAL(n))*(hi - lo);
    table(i,2) = P(table(i,1));
    table(i,3) = exp(table(i,1));
    table(i,4) = log(table(i,2));
  }

  writetable("yPCp.txt",table,header,15,10);

  header[0] = "y";
  header[1] = "dp";
  header[2] = "p";
  header[3] = "dq";

  for (INTEGER i=1; i<=n; ++i) {
    table(i,1) = lo + (REAL(i)/REAL(n))*(hi - lo);
    table(i,2) = p.derivative(table(i,1));
    table(i,3) = p(table(i,1));
    table(i,4) = q.derivative(table(i,3));
  }

  writetable("dPdQ.txt",table,header,15,10);

  n = 1000;

  INT_32BIT seed = 740726;

  REAL sdev_y = sigma/sqrt(1.0 - alpha*alpha);
  REAL y_lag = sdev_y*unsk(seed);

  realmat pc(n,2);

  for (INTEGER i=1; i<=n; ++i) { 
    REAL y = alpha*y_lag + sigma*unsk(seed);
    pc(i,1) = p(y);
    pc(i,2) = y;
    y_lag = y;
  }

  writetable("pc.txt",pc,20,16);

  return 0;
}
