#!/usr/bin/python
"""
This program will simulate bond prices for the cash flow matching problem.

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 3 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, see <http://www.gnu.org/licenses/>.
"""

import sys
import math
import random
from optparse import OptionParser

def P0(T, param):
    return math.exp(param.B / param.C * (math.exp(- param.C * T)-1)  - param.A * T)

def A(t, T, param):
    return math.log(P0(T, param) / P0(t, param)) \
            + B(t, T, param) * F(t, param) - \
            param.sigma ** 2 / 4 / param.alpha * \
            (1 - math.exp(- param.alpha * (T - t))) ** 2 * \
            (1 - math.exp(- 2 * param.alpha * t))

def B(t, T, param):
    return (1 - math.exp(- param.alpha * (T - t))) / param.alpha

def F(t, param):
    return param.A + param.B * math.exp(- param.C * t)

def g(u, t, param):
    return F(t, param) - F(u, param) * math.exp(-param.alpha * (t - u)) \
            + param.sigma ** 2 / 2 / param.alpha * ( 
                (1 - math.exp(-param.alpha * t)) ** 2 - 
                (1 - math.exp(-param.alpha * u)) ** 2 * 
                math.exp(-param.alpha * (t - u))
                )

def main():
    usage = "usage: %prog [option] coupon file"
    version = "%prog 1.0"
    parser = OptionParser(usage=usage, version=version)
    parser.add_option("--horizon", "-T", dest="T", action="store",
            help="time horizon (T) in year", default=0, type = 'int')
    parser.add_option("--frequency", "-q", dest="freq", action="store",
            help="number of time steps in a year, (>= 1)", 
            default=4, type = 'int')
    parser.add_option("--alpha", "-a", dest="alpha", action="store",
            help="alpha", default=None, type = 'string')
    parser.add_option("--sigma", "-s", dest="sigma", action="store",
            help="sigma", default=None, type = 'string')
    parser.add_option("--constantA", "-A", dest="A", action="store",
            help="A as in F(t) = A + B e^{-C t}", 
            default=None, type = 'string')
    parser.add_option("--constantB", "-B", dest="B", action="store",
            help="B as in F(t) = A + B e^{-C t}", 
            default=None, type = 'string')
    parser.add_option("--constantC", "-C", dest="C", action="store",
            help="C as in F(t) = A + B e^{-C t}", 
            default=None, type = 'string')
    parser.add_option("--samplesize", "-N", dest="N", action="store",
            help="Number of samples", default=1, type = 'int')
    (option, arg) = parser.parse_args()

    # check the options
    try:
        option.alpha = float(option.alpha)
    except:
        parser.error("please check the alpha (-a)")
    try:
        option.sigma = float(option.sigma)
    except:
        parser.error("please check the sigma (-s)")
    try:
        option.A = float(option.A)
    except:
        parser.error("please check the A (-A)")
    try:
        option.B = float(option.B)
    except:
        parser.error("please check the B (-B)")
    try:
        option.C = float(option.C)
    except:
        parser.error("please check the C (-C)")

    # get the coupon
    coupon = []
    if arg == []:
        for i, j in enumerate(sys.stdin):
            coupon.append(map(float, j.split()))
    elif len(arg) == 1:
        f = open(arg[0], 'r')
        for i, j in enumerate(f):
            coupon.append(map(float, j.split()))
        f.close()
    else:
        return 0


    # simulate r(t)
    # r(0), r(1), ..., r(N - 1)
    sim = 0
    delta_t = 1.0 / option.freq 
    while sim < option.N:
        r = []
        r.append(F(0, option))
        for i in range(1, option.T * option.freq):
            r.append(
                    math.exp(- option.alpha * delta_t) * r[-1] + 
                    g((i - 1 ) * delta_t, i * delta_t, option) +
                    math.sqrt(option.sigma ** 2 / 2 / option.alpha * (
                        1 - math.exp(-2 * option.alpha * delta_t)
                        )) * random.normalvariate(0, 1)
                    )
        for t, i in enumerate(r):
            for j in coupon:
                price = 0.0
                for k, l in enumerate(j):
                    if l != 0:
                        price += l * math.exp(
                                A(t * delta_t, (t + k + 1) * delta_t, option) - 
                                B(t * delta_t, (t + k + 1) * delta_t, option) 
                                * i)
                print price,
        sim += 1
        if (sim != option.N):
            print '\n',

if __name__ == "__main__":
    main()