/*
 * diff.cpp
 *
 * Created on April 11, 2006, 2:25 PM
 */
#ifndef DIFF_CPP
#define DIFF_CPP
/**
 *
 * @author horacio a. aliaga
 */
#include <vector>
#include <iostream>
#include <fstream>
#include <string>
#include <math.h>
#include "SdeIntegrator.h"
#include "TwoCurves.cpp"

using namespace std;
/** BDT interest rate model in differential form, where  
 * r is the short term risk free rate of return, 
 * dt is the time step differential, 
 * r is pulled back to the mu-level at a rate of alpha, 
 * a represents the function a(r)=alpha(mu-sigma'/sigma ln r) used in 
 * the SDEIntegrator. The function b(r) is the time dependent volatility. 
 * The sde is an SDEIntegrator used to find r(t+delta_t).
 */
class diff {
public: 
	double r, r0;
    double B_maturity; //Bond maturity 
    double dt;
    double a, d_a,dd_a;
    double b, d_b, dd_b, ddd_b;
    double* time;
    double* yield;
    double* yield_prime;
    double* sigma;
    double* sigma_prime;
    int index;
    double call,put;
    SDEIntegrator sde;
    
    /**
     * Creates a new instance of diff with 
     * default values
     */
public: diff() {
        dt=0.0;
        r = 0.0;
        r0 = 0.0;
        index =0;
    }
    /**
     * Creates a new instance of diff
     */
public: diff(double _dt,double* _time,double* _yield, 
			 double* _sigma,double* _sigma_prime, double _r0) {
        dt=_dt;
        time=_time;
        yield=_yield;
        sigma=_sigma;
        sigma_prime=_sigma_prime;
        r = _r0;
        r0 = _r0;
        index =0;
        set_b();set_d_b();set_dd_b();set_ddd_b();
        set_a();set_d_a();set_dd_a();
        sde = SDEIntegrator(a,d_a,dd_a,b,d_b,dd_b,ddd_b,dt,r);
    }
public: void set_SDE() {
        set_b();set_d_b();set_dd_b();set_ddd_b();
        set_a();set_d_a();set_dd_a();
        sde.set_SDEIntegrator(a,d_a,dd_a,b,d_b,dd_b,ddd_b,dt,r);
    }
public: void set_B_maturity(double b){
        B_maturity=b;
    }
public: double get_B_maturity(){
        return B_maturity;
    }
public: void set_index(int i){
        index=i;
    }
public: int get_index(){
        return index;
    } 
    /** sets the a(r) function used in the SDEIntegrator*/
public: virtual void set_a(){};
public: virtual void set_d_a(){};
public: virtual void set_dd_a(){};
    /** sets the b(r,t) function used in the SDEIntegrator*/
public: virtual void set_b(){};
    /** sets the db(r,t)/dr function used in the SDEIntegrator*/
public: virtual void set_d_b(){};
    /** sets the d^2b(r,t)/dr^2 function used in the SDEIntegrator*/
public: virtual void set_dd_b(){};
    /** sets the d^3b(r,t)/dr^3 function used in the SDEIntegrator*/
public: virtual void set_ddd_b(){};
public: void set_call(double c){
        call=c;
    }
public: void set_put(double p){
        put=p;
    }
public: void set_dt(double _dt){
        dt=_dt;
    }
public: double get_dt(){
        return dt;
    }
    /** gets a(r)*/
public: double get_a(){
        return a;
    }
public: double get_d_a(){
        return d_a;
    }
public: double get_dd_a(){
        return dd_a;
    }
    /** sets the short term risk free interest rate */
public: void set_r(double _p){
        r=_p;
    }
    /** sets the initial short term risk free interest rate */
public: void set_r0(double _p){
        r0=_p;
    }
    /** sets the rate at which r(t) goes back to mu*/
public: void set_time(double* _time){
        time=_time;
    }
    /** long term average interest rate */
public: void set_yield(double* _yield){
        yield=_yield;
    }
    /** long term average interest rate */
public: void set_yield_prime(double* _yield_p){
        yield_prime=_yield_p;
    }
    /** sets the volatility sigma of r(t)*/
public: void set_sigma(double* _sigma){
        sigma=_sigma;
    }
    /** sets the first derivative of the volatility sigma respect of t*/
public: void set_sigma_prime(double* _sigma_prime){
        sigma_prime=_sigma_prime;
    }
    /** gets the short term risk free interest rate*/
public: double get_r(){
        return r;
    }
    /** gets the initial short term risk free interest rate*/
public: double get_r0(){
        return r0;
    }
    /** gets the rate alpha at which r(t) converges to mu*/
public: double* get_time(){
        return time;
    }
    /** gets the long term average rate of return mu*/
public: double* get_yield(){
        return yield;
    }
    /** gets the long term average rate of return mu*/
public: double* get_yield_prime(){
        return yield_prime;
    }
    /** gets the volatility sigma of r(t)*/
public: double* get_sigma(){
        return sigma;
    }
    /** gets the d(sigma)/dt  */
public: double* get_sigma_prime(){
        return sigma_prime;
    }
public: double get_b(){
        return b;
    }
public: double get_d_b(){
        return d_b;
    }
public: double get_dd_b(){
        return dd_b;
    }    
public: double get_ddd_b(){
        return ddd_b;
    }    
    /** The Polar Marsaglia generator of rand gaussian variables*/
private: double pol_mas(double t)
    {
	double u1,u2,v1,v2,w,wn,dwt,dwt2;
        dwt =0.0;
        for(int i=0;i<10000;i++)
        {
            u1=rand();
            v1=2*u1-1;
            u2=rand();
            v2=2*u2-1;
            w=v1*v1+v2*v2;
            if ((w <= 1.0)&&(w >0.0))
            {
                wn=sqrt(-2*log(w)/w);	  
                dwt=sqrt(t)*v1*wn;
                dwt2=sqrt(t)*v2*wn;
                break;
            }
        }
        return dwt;
    }    
   /** Curve[0] has the time dependent IR
     * Curve[1] has the time dependent price of the Bond
     */    
public: TwoCurve gen_BondandIRvsMaturity(int n_points,double S){
        // TODO code application logic here
        TwoCurve curve;
        double y,y0,y1,y2;
        double* r=new double[n_points+1];
        double* bond=new double[n_points+1];

        double t=0;
        sde.set_X(sde.get_X0());
        set_r(get_r0());
        y = log(get_r0());
        double r_aver =get_r();
        set_index(0);
        for (int i=0;i<n_points;i++){        
            r[i]=r[i]+exp(y);
            double ax=exp(-r_aver/(1.0*(i+1))*t);
            bond[i]=bond[i]+S*ax;
            Point* a = new Point(t,get_r());
            curve.c0_add(a);
            Point* b = new Point(t, S*ax);
            curve.c1_add(b);
            y0=sde.Strong_Integrator_Zero();
            y1=sde.Strong_Integrator_First();
            y2=sde.Strong_Integrator_Second();
            //y3=v.sde->Strong_Integrator_Third();
            //y4=v.sde->Strong_Integrator_Fourth();
            y=y0+y1+y2;
            set_r(exp(y));
            set_index(i+1);
            set_b();set_d_b();set_dd_b();set_ddd_b();
            set_a();set_d_a();set_dd_a();
            set_SDE();sde.set_X(exp(y));
            r_aver = r_aver + exp(y);
            t = t + get_dt();
        }
        y0=sde.Strong_Integrator_Zero();
        y1=sde.Strong_Integrator_First();
        y2=sde.Strong_Integrator_Second();
        y=y0+y1+y2;
        r_aver = r_aver + exp(y);;
        r[n_points]=r[n_points]+exp(y);
        double ax=exp(-r_aver/(1.0*(n_points))*t);
        bond[n_points]=bond[n_points]+S*ax;
        return curve;
    }            
public: TwoCurve gen_AverBondandIRvsMaturity(int n_sim,int n_points,double S) {
        // TODO code application logic here
        double y,y0,y1,y2;
        double* r = new double[n_points+1];
        double* bond = new double[n_points+1];
        for (int i=0;i<n_points;i++){
			r[i]= 0.0;
			bond[i]=0.0;
		}
        for(int j=0;j<n_sim;j++){
            double t=0;
            sde.set_X(sde.get_X0());
            set_r(get_r0());
            y = log(get_r0());
            double r_aver =get_r();
            set_index(0);
            for (i=0;i<n_points;i++){
                r[i]=r[i]+exp(y);
                double ax=exp(-r_aver/(1.0*(i+1))*t);
                bond[i]=bond[i]+S*ax;
                y0=sde.Strong_Integrator_Zero();
                y1=sde.Strong_Integrator_First();
                y2=sde.Strong_Integrator_Second();
                //y3=v.sde->Strong_Integrator_Third();
                //y4=v.sde->Strong_Integrator_Fourth();
                y=y0+y1+y2;
                set_r(exp(y));
                set_index(i+1);
                set_b();set_d_b();set_dd_b();set_ddd_b();
                set_a();set_d_a();set_dd_a();
                set_SDE();sde.set_X(exp(y));
                r_aver = r_aver + exp(y);
                t = t + get_dt();
            }
            y0=sde.Strong_Integrator_Zero();
            y1=sde.Strong_Integrator_First();
            y2=sde.Strong_Integrator_Second();
            y=y0+y1+y2;
            r_aver = r_aver + exp(y);;
            r[n_points]=r[n_points]+exp(y);
            double ax=exp(-r_aver/(1.0*(n_points))*t);
            bond[n_points]=bond[n_points]+S*ax;
        }
        TwoCurve curve;
        double t=0;
        for (i=0;i<=n_points;i++){
            r[i]=r[i]/(1.0*n_sim);
            bond[i]=bond[i]/(1.0*n_sim);
            Point* a = new Point(t,r[i]);
            curve.c0_add(a);
            Point* b = new Point(t, bond[i]);
            curve.c1_add(b);
            t = t + get_dt();
        }
        return curve;
    }                
};

#endif DIFF_CPP