##
## The data: 
##
BQ = data.frame(coupon=c(6 + 5/8, 9 + 1/8, 7 + 7/8, 8 + 1/4, 8 + 1/4, 8 + 3/8, 8, 8 + 3/4, 6 + 7/8, 8 + 7/8, 6 + 7/8, 8 + 5/8, 7 + 3/4, 11 + 1/4, 8 + 1/2, 10 + 1/2, 7 + 7/8, 8 + 7/8),
                ask_price=c(100 + 0/32, 100 + 22/32, 100 + 24/32, 101 + 1/32, 101 + 7/32, 101 + 12/32, 100 + 26/32, 102 + 1/32, 98 + 5/32, 102 + 9/32, 97 + 13/32, 101 + 23/32, 99 + 5/32, 109 + 4/32, 101 + 13/32, 107 + 27/32, 99 + 13/32, 103 + 0/32))

## It looks like the semiannual coupon payments will be paid on:
##
## February 15th
## August 15th
##
## each year.

## Calculate the accrued interest between the-last coupon payment (on 20110815) and today (on 20111105)
##
## As there are 92 days between these two dates and 92 days until the next coupon payment the current
## owner of the bond is owed this fraction of the half-coupon:
##
## See Page 51 of the book "Investment Science" on calculating accrued interest.
##
ndays_from = 82
ndays_till = 102

accrued_interest = (ndays_from/(ndays_from + ndays_till)) * (BQ$coupon/2)

BQ$ask_price_w_accrued_interest = BQ$ask_price + accrued_interest

## What numbers did the book compute:
BQ$books_ask_price_w_accrued_interest = c(101.48, 102.72, 102.5, 102.87, 103.06, 103.24, 102.59, 103.99, 99.69, 104.26, 98.94, 103.63, 100.9, 111.64, 103.3, 110.18, 101.17, 104.97)

if( FALSE ){
    ## How well do my numbers match the book's:
    plot (BQ$books_ask_price_w_accrued_interest, BQ$ask_price_w_accrued_interest, type='p')
    abline(a=0, b=1, col='green')
    grid()
}

## The first coupon paid to us is smaller (since some should have already been payed out to the previous owner of the bond):
##
frac_first_coupon = ndays_till/(ndays_from + ndays_till) ## fraction of the first coupon we will get
x0 = frac_first_coupon * BQ$coupon/2 ## the amount to be recieved on the first coupon payment on 20111215
t0 = ndays_till/365. ## time in years until first coupon on 20111215

## The time in years until expiration:
##
BQ$time_to_maturity = c(rep(102, 2), rep(284, 2), rep(468, 2), rep(649, 2), rep(833, 2), rep(1014, 2), rep(1198, 2), rep(1379, 2), rep(1563, 2))/365

discount_rate = function(t, a0, a1, a2, a3, a4){
    a0 + a1*t + a2*t*2 + a3*t^3
}

## Compute the present values of all of the coupons payed until maturity using a discount rate function:
##

compute_present_values = function(BQ, t0, x0, a0, a1, a2, a3, a4){
    ## Process the first coupon date:
    ##
    t = t0
    x = x0
    PV = x0 * exp(-t*discount_rate(t, a0, a1, a2, a3))

    bonds_ending = (BQ$time_to_maturity > t - 0.25) & (BQ$time_to_maturity < t + 0.25)
    ##print('bonds_ending')
    ##print(bonds_ending)
    PV = PV + 100 * bonds_ending * exp (-t*discount_rate(t, a0, a1, a2, a3, a4))

    bonds_ended = BQ$time_to_maturity < t + 0.25
    ##print('bonds_ended')
    ##print(bonds_ended)
    
    while( any(!bonds_ended) ){
        t=t+0.5 ## the time of the next coupon i.e. 6 months

        bonds_paying_coupons = (BQ$time_to_maturity > t - 0.25)
        ##print('bonds_paying_coupons')
        ##print(bonds_paying_coupons)
        PV = PV + bonds_paying_coupons * ( BQ$coupon/2 ) * exp(-t*discount_rate(t, a0, a1, a2, a3, a4))

        bonds_ending = (BQ$time_to_maturity > t - 0.25) & (BQ$time_to_maturity < t + 0.25)
        ##print('bonds_ending')
        ##print(bonds_ending)
        PV = PV + 100 * bonds_ending * exp (-t*discount_rate(t, a0, a1, a2, a3, a4))

        bonds_ended = BQ$time_to_maturity <t + 0.25
        ##print('bonds_ended')
        ##print(bonds_ended)
    }
    PV
}

## Test the above function (using what I think are parameters close to the optional ones):
##
print(compute_present_values(BQ, t0, x0, a0=0.062009143, a1=0.00627032, a2=0.001099467, a3=-0.000593607, a4=4.95e-6))

## Here is a function that one could use to implement an optimization routine using:
##
min_fn = function(a0, a1, a2, a3, a4, BQ=BQ, t0=t0, x0=x0){
    PV_hat = compute_present_values(BQ, t0, x0, a0, a1, a2, a3, a4)
    sum((BOQ$ask_price_w_accrued_interest - PV_hat)^2)
}