/* \file generalizedOUprocess.hpp
   \brief Ornstein-Uhlenbeck process with piece-wise linear coefficients
*/
#include "..\generalizedOUprocess.hpp"

namespace QuantLib {

	GeneralizedOUprocess::GeneralizedOUprocess(const Handle<YieldTermStructure>& speedTS,
		const Handle<YieldTermStructure>& vol,
		Real x0,
		Real level)
		: x0_(x0), speed_(speedTS), level_(level), volatility_(vol) {

					QL_REQUIRE(x0 >= 0.0, "negative initial data given");
			QL_REQUIRE(level >= 0.0, "negative level given"); 
	}

	Real GeneralizedOUprocess::x0() const {
		return x0_;
	}

	Handle<YieldTermStructure> GeneralizedOUprocess::speed() const {
		return speed_;
	}

	Handle<YieldTermStructure> GeneralizedOUprocess::volatility() const {
		return volatility_;
	}
	
	Real GeneralizedOUprocess::drift(Time t, Real x) const {
        
		Real speed;
			if (t!=0) {
				speed = - std::log(speed_->discount(t))/t;
			} else { speed = 1.0;}

        return speed * (level_ - x);;
    }

	Real GeneralizedOUprocess::diffusion(Time t, Real) const {

        Volatility vol;
		if (t!=0) {
				vol = - std::log(volatility_->discount(t))/(t);
			} else {
				vol = 0.3;
		}

        return vol;
    }

	Real GeneralizedOUprocess::level() const {
		return level_;
	}

	Real GeneralizedOUprocess::expectation(const Time t, Real x0, Time dt) const {

		Real speed;
			if (t!=0) {
				speed = - std::log(speed_->discount(t))/t;
			} else {
				// to avoid division by zero
				speed = 1.0;
			}

			return level_ + (x0 - level_) * std::exp(-speed*dt);
	}

	Real GeneralizedOUprocess::stdDeviation(Time t, Real x0, Time dt) const {
			return std::sqrt(variance(t,x0,dt));
	}

	Real GeneralizedOUprocess::variance(const Time t, Real x, Time dt) const {
			Real speed;
			Volatility vol;
		if (t!=0) {
				speed = - std::log(speed_->discount(t))/t;
				vol =  - std::log(volatility_->discount(t))/(t);
			} else {
				// To make sure that the tree opens up nicely at the start and to avoid division by zero
				speed = 1.0;
				vol = 0.3;
		}

		if (speed < std::sqrt(QL_EPSILON)) {
			// algebraic limit for small speed
			return vol*vol*dt;
		} else {
			return 0.5*vol*vol/speed*
				(1.0 - std::exp(-2.0*speed*dt));
		}
	}

}

