423.sagews

# C : periodic coupon (interest rate) [ex, 5% would be .05]
# n : number of periods [ex, 10]
# F : maturity value [ex, 1000]
# r : required yield [ex, 5% would be .05]
def compute_units(C, n, F, r):
    return C * F * ((1 - (1 / (1 + r) ^ n)) / r) + F * (1 / (1 + r) ^ n)

# Same as compute_units, except returns a scaled price (100 being par)
# instead of the units.
def compute_price(C, n, F, r):
    return 100 * (compute_units(C, n, F, r) / F)

# Graph a price yield curve.  7% coupon, $1000 principal, with yield ranging from 0% to 20%.
# The red curve is a 10 year bond, and the blue is a 1 year bond.
p = plot(compute_price(.07,1,1000,var('r')),0,.2)
p += plot(compute_price(.07,10,1000,var('r')),0,.2, rgbcolor=(1,0,0))
show(p)

# Same scenario as before, but showing the time periods as a 3rd dimension instead of
# overlaying the curves.  Note: the previous 2 curves are the endpoints of this surface.
plot3d(compute_price(.07,var('n'),1000,var('r')), (n,1,10), (r,0,.2))

# Newton-Raphson function.  Takes a function, its derivative, an initial guess, and 
# an optional tolerance.
def newton(func, funcd, x, TOL=1e-6):
    f, fd = func(x), funcd(x)
    count = 0
    while 1:
        dx = f / float(fd)
        if abs(dx) < TOL * (1 + abs(x)): return x - dx
        x = x - dx
        f, fd = func(x), funcd(x)
        count = count + 1
        print "newton(%d): x=%s, f(x)=%s" % (count, x, f)

# Compute the yield given the bond characteristics, a price, and an initial guess.
def compute_yield(C, n, F, u, y0):
    r = var('r')
    func = compute_units(C,n,F,r) - u
    funcd = diff(func, r)
    return newton(func, funcd, y0)

C = .0675   # coupon
n = 20      # periods
F = 10000   # principal
u = 9125    # bond price (units)
g = .1      # initial guess for yield

# This is like computing f(g(f(x))), where f computes units and g computes yield.
compute_units(C, n, F, compute_yield(C, n, F, u, .1))

newton(1): x= .07144291315814208, f(x)= 461.9310522332926 newton(2): x= .07599073225081789, f(x)= 15.89642112795082 newton(3): x= .07615863583051595, f(x)= .02048624306917191 9125.00000003