Yield To Maturity calculation using a Python script

When interested in buying bonds a term that often comes up is the definition of Yield To Maturity (YTM):

Yield to maturity (YTM) is the total return anticipated on a bond if the bond is held until the end of its lifetime. Yield to maturity is considered a long-term bond yield, but is expressed as an annual rate. In other words, it is the internal rate of return of an investment in a bond if the investor holds the bond until maturity and if all payments are made as scheduled. – Investopedia

Let us consider a bond with a par value of \$100 and a coupon rate of 3\% with the bond maturing in 29 years time. The current bond price is \$89. We are interested in the bond’s YTM.

Knowing that we get \$3 annual cashflow (coupon\ rate \times par\ value = 0.03 \times \$100 = \$3) we will have the following cashflow for the next 29 years. Notice in the last year we will receive an additional \$100 which is the bond’s principle:

Let us assume that we want to value these future cashflows and determine how much they are worth today. What we need to do, is to present value all of these future cashflows back to today and add up those present values pv_i. This sum would be the price of the bond.

A present value is the amount of money pv one would need to invest today in order to have c amount of money in n years if he is earning r% (annual)interest on his investment:

  pv = \frac{c}{(1 + r)^n}

For example, if Bob wants to have \$100 one year from now at a rate of 5\%, he needs to invest \$95.24 today (\frac{100}{(1 + 0.05)} = 95.24).

Back to our bond’s present value of future cashflows the resulting formula would look as follows:
  bond\ price = \frac{future\ cashflow\ 1}{(1 + ytm)^1} + \frac{future\ cashflow\ 2}{(1 + ytm)^2} + \frac{future\ cashflow\ 3}{(1 + ytm)^3} + \dots + \frac{future\ cashflow\ 29}{(1 + ytm)^{29}}

with ytm being the yield to maturity we are looking for.

For our particular bond the formula will look like this:
  bond\ price = \frac{3}{(1 + ytm)^1} + \frac{3}{(1 + ytm)^2} + \frac{3}{(1 + ytm)^3} + \dots + \frac{3 + 100}{(1 + ytm)^{29}}

So how can we go about finding ytm? Well, above we have learned that sum of present values of all future cash flows is our bond price. Since we already know the today’s bond price which is \$89, the formula becomes like this:
  \frac{3}{(1 + ytm)^1} + \frac{3}{(1 + ytm)^2} + \frac{3}{(1 + ytm)^3} + \dots + \frac{3 + 100}{(1 + ytm)^{29}} = 89

All we need to do now is to solve the above equation for ytm so that it holds true. We will do this through trial-and-error.

Solving the equation by hand requires an understanding of the relationship between a bond’s price and its yield, as well as of the different types of bond pricings. Bonds can be priced at a discount, at par and at a premium. When the bond is priced at par, the bond’s interest rate is equal to its coupon rate. A bond priced above par (called a premium bond) has a coupon rate higher than the interest rate, and a bond priced below par (called a discount bond) has a coupon rate lower than the interest rate. So if an investor were calculating YTM on a bond priced below par, he or she would solve the equation by plugging in various annual interest rates that were higher than the coupon rate until finding a bond price close to the price of the bond in question (quote from: Investopedia).

You can verify this quickly with a simple example. Assume a bond with current price b and a par value of 100 that matures in 1 year. The coupon rate is c and interest rate is ytm:

  b = \frac{100 + c \times 100}{(1 + ytm)}
==>
  b \times (1 + ytm) = 100 + c \times 100 = 100 \times (1 + c)
==>
  \frac{b}{100} = \frac{(1 + c)}{(1 + ytm)}

If \frac{b}{100} \textless 1 then (1 + ytm) \textgreater (1 + c). Similarly for the other cases.

Back to our bond. We had the following equation:
  \frac{3}{(1 + ytm)^1} + \frac{3}{(1 + ytm)^2} + \frac{3}{(1 + ytm)^3} + \dots + \frac{3 + 100}{(1 + ytm)^{29}} = 89

Now we must solve for the interest rate ytm which is where things start to get difficult. Yet, we do not have to start simply guessing random numbers if we stop for a moment to consider the relationship between bond price and yield. As was mentioned above, when a bond is priced at a discount from par, its interest rate will be greater than the coupon rate. In this example, the par value of the bond is \$100, but it is priced below the par value at \$89, meaning that the bond is priced at a discount. As such, the annual interest rate we are seeking must necessarily be greater than the coupon rate of 3\%.

With this information we can now calculate and test a number of bond prices by plugging various annual interest rates that are higher than 3\% into the formula above. Using a few different interest rates above 3\%, one would come up with the following bond prices:

I have used 0.00001 increments and as we can see with an interest rate of 3.62\% we can approximate our current price of \$89 quite good enough. So 3.62\% would be our Yield To Maturity.

Full code

References
Yield To Maturity (YTM) – Investopedia

 

The Knuth shuffle

The Knuth shuffle algorithm is an algorithm for randomly shuffling the elements of an array arr of size n given a function that generates a uniformly distributed random number in \mathcal{O}(1) time. The algorithm produces an unbiased permutation. That is, every permutation is equally likely.

Proof that the algorithm is unbiased
For any element e, the probability that it will be shuffled into position n-1 (last position at zero based indexing) is: the probability that it gets selected when i=n-1, which is \frac{1}{n}.

For any element e, the probability that it will be shuffled into position n-2 (second-to-last position) when i=n-2 is: the probability that it gets selected when i=n-2, which is \frac{1}{n-1} multiplied by the probability that it is not selected for the last position, which is 1 - \frac{1}{n} = \frac{n-1}{n}, giving  \frac{1}{n-1}\times\frac{n-1}{n} = \frac{1}{n}

For any element e, the probability that it will be shuffled into position n-3 (third-to-last position) when i=n-3 is: the probability that it gets selected when i=n-3, which is \frac{1}{n-2} multiplied by the probability that it is not selected for the last nor second-to-last position, which is 1 - \frac{2}{n} = \frac{n-2}{n}, giving \frac{1}{n-2}\times\frac{n-2}{n} = \frac{1}{n}.

This can easily be generalised for any other position.

So we can see that for any element e, the probability that it will be shuffled into the last, second-to-last, third-to-last, etc position is \frac{1}{n}.

 

Implement a ceil() function

In this post I will explain a less intuitive approach of implementing a ceil() function in Java, that is faster than known comparison approaches.

A few known facts
Fact 1
Dividing two integers a and b in Java will always floor the result:

well almost. Actually Math.floor(a/b) rounds towards negative infinity, where as a/b rounds towards zero:

Fact 2
\lfloor x \rfloor = m if m \le x \textless m+1
\lceil x \rceil = n if n - 1 \textless x \le n

Fact 3
\lfloor x \rfloor \le \lceil x \rceil

What we need is \lfloor (x + s) \rfloor == \lceil x \rceil
Since dividing two integers a and b in Java will always floor the result, we need to find some s such that \lfloor (x + s) \rfloor == \lceil x \rceil.

How do we find s?
Think about it this way: \lceil \frac{a}{b} \rceil \textgreater \lfloor \frac{a}{b} \rfloor, except when \frac{a}{b} is a whole number. So, we want to bump \frac{a}{b} to (or past) the next whole number, unless \frac{a}{b} is a whole number, by adding s. This way \lfloor \frac{a}{b} + s \rfloor == \lceil \frac{a}{b} \rceil. We want to find s so that for whole numbers, it won’t quite push them up to the next whole number, but for non-whole numbers it will.

For instance, assume a=6, b=3 and s=1 then \lfloor \frac{a}{b} + s \rfloor = \lfloor \frac{6}{3} + 1 \rfloor = 3 \ne \lceil \frac{a}{b} \rceil = \lceil \frac{6}{3} \rceil = 2, which is not the desired result. As you can see, by using s=1 we bumped \frac{a}{b} to far. So s must definitely be less than 1.

So how much s is just enough?

Assume a=7, b=3. Then \frac{a}{b} = \frac{7}{3}, which can be re-written as:
\frac{7}{3} = \frac{1}{b} + \frac{1}{b} + \frac{1}{b} + \frac{1}{b} + \frac{1}{b} + \frac{1}{b} + \frac{1}{b} = (\frac{1}{3} + \frac{1}{3} + \frac{1}{3}) + (\frac{1}{3} + \frac{1}{3} + \frac{1}{3}) + \frac{1}{3} = 2  + \frac{1}{3}
In other words, we get an integer part (2) and a left over (\frac{1}{3}). The left over part is always in the form of k\frac{1}{b}, with 0 \le k \textless b. Thus, the smallest non-whole number we can get from \frac{a}{b} is qb + \frac{1}{b} with qb being the integer part with q \in \mathbb{N}_0.

So assuming \frac{a}{b} is not a whole number, the smallest left over we can have is \frac{1}{b}. So in-order to bump \frac{a}{b} to the next whole number, it suffices to add s = 1 - \frac{1}{b}. This is \textless 1, which will not bump \frac{a}{b} to the next whole number in case \frac{a}{b} is a whole number, but still enough to bump \frac{a}{b} to the next whole number in case \frac{a}{b} is not a whole number.

Summing it up
So we know adding s = 1 - \frac{1}{b} to \frac{a}{b} will fulfill the equality: \lfloor (\frac{a}{b} + s) \rfloor = \lceil \frac{a}{b} \rceil. Thus we get:

 \lceil \frac{a}{b} \rceil = \lfloor (\frac{a}{b} + s) \rfloor = \lfloor (\frac{a}{b} + 1 - \frac{1}{b}) \rfloor = \lfloor \frac{a + b - 1}{b} \rfloor

 

Exponentiation by squaring

This is a rather short post. in which I explain how to we can implement a fast pow function to raise a number x to the power of n.

A rather naive approach would be to do the following:

This is ok, but the algorithm has a running time of \mathcal{O}(n). We can do better than that.

Let us consider an example.

Assume that we want to calculate x^8. We can rewrite it as follows:

x^8 = x \times x \times x \times x \times x \times x \times x \times x.

which requires 7 multiplications. What we can do to reduce the number of multiplications is to first compute:

x \times x = x^2

next we can compute

x^2 \times x^2 = x^4

and last

x^4 \times x^4 = x^8

which gives as a total of 3 multiplications.

Lets assume that we want to calculate x^{13}. We can rewrite it as follows:

x^{13} = x^8 \times x^4 \times x.

with 8, 4 and 1 being all powers of 2. The great thing about 8, 4 and 1 being a power of 2 is that we can easily get them by squaring (see previous example).

How do we find the numbers 8, 4 and 1? Or in other words, how can we find out what powers of 2 make up 13?
Well that’s what the binary representation of 13 tells us:

13_d = 1101_b

Now that we have our powers of 2, all we need to do is to keep squaring (that is compute x \times x, then (x \times x) \times (x \times x) , etc) and every time we reach a power of two (‘1’ in binary representation) we multiply the intermediate results. So in order to compute x^{13} we iterate over the binary representation of 13 from the least significant bit (right-most bit) to left-most and perform our squaring and multiplication.

Algorithm

  1. Set result = 1
  2. 1101  \rightarrow we have a power of 2 (2^0), so we multiply result with x: result = result \times x
  3. Square x: x = x \times x  \rightarrow x = x^2
  4. 1101  \rightarrow we don’t have a power of 2 so skip.
  5. Square x: x = x \times x  \rightarrow x = x^4
  6. 1101  \rightarrow we have a power of 2 (2^2), so we multiply result with x: result = result \times x
  7. 1101  \rightarrow we have a power of 2 (2^3), so we multiply result with x: result = result \times x

Code

The running time thus becomes \mathcal{O}(\log n)