James Gray
2/11/2008 12:25:00 AM
[Note: parts of this message were removed to make it a legal post.]
Begin forwarded message:
> From: "Clark Grubb" <clarkgrubb@gmail.com>
> Date: February 10, 2008 1:45:10 AM CST
> To: submission@rubyquiz.com
> Subject: Please Forward: Ruby Quiz Submission
>
> Since it's possible for a cash flow to have multiple IRRs, this code
> returns an array
> containing all the solutions. We do this by finding the roots of
> the polynomial in (1+IRR)
> obtained by clearing the IRR equation of denominators. We discard
> roots with imaginary component
> or where (1+IRR) <= 0.
>
> To find the roots, we use Laguerre's method, which is described in
> "Numerical Recipes in C++",
> 2nd edition, pp. 376-379. Although Laguerre's is more complex to
> implement than Newton's method
> or bisection, it is easier to use since it converges reliably
> regardless of starting point.
>
> =========================
>
> Require 'Complex'
>
> class Polynomial
>
> attr_accessor :epsilon, :max_iterations
>
> @@laguerre_fractions = [0.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62,
> 0.88, 1.0]
>
> def initialize(coefficients)
> @coefficients = coefficients
> @epsilon = 1e-12
> @max_iterations = 80
> end
>
> def [](i)
> @coefficients[i] || 0
> end
>
> def []=(i,v)
> @coefficients[i]=v
> end
>
> def degree
> d = 0
> @coefficients.each_with_index { |c,i| d = i unless c.nil? or
> c.zero? }
> @coefficients = @coefficients.slice(0..d)
> @coefficients.size - 1
> end
>
> def deflate!(z)
> b = 0
> degree.downto(0) { |i| self[i],b = b, z*b + self[i] }
> end
>
> def laguerre(z=nil)
> z ||= Complex.new(0.0,0.0)
> i_per_frac = @max_iterations / (@@laguerre_fractions.size - 1)
> (i_per_frac * (@@laguerre_fractions.size - 1)).times do |i|
> b = self[degree]
> err = b.abs
> d = f = 0.0
> m = degree
> (m-1).downto(0) do |j|
> f = z*f + d
> d = z*d + b
> b = z*b + self[j]
> err = b.abs + z.abs * err
> end
> err *= epsilon
> return z if b.abs <= err
> g = d/b
> h = g*g - 2.0 * f / b
> sq = Math.sqrt( (m-1)*(m*h -g*g) )
> gp = g+sq
> gm = g-sq
> gp = gm if gp.abs < gm.abs
> if [gp.abs,gm.abs].max > 0.0
> dz = m/gp
> else
> dz = (1+z.abs) * Math::exp( Complex::I * i )
> end
> return z if z == z-dz
> if 0 != i % i_per_frac
> z = z-dz
> else
> z -= @@laguerre_fractions[i/i_per_frac]*dz
> end
> end
> raise "failed to converge after #{max_i} iterations"
> end
>
> def roots
> if @roots.nil?
> @roots = []
> p = self.dup
> degree.times do
> z = p.laguerre
> if z.image.abs <= 2 * epsilon * z.real.abs
> z = Complex.new(z.real,0)
> end
> @roots << z
> p.deflate!(z)
> end
> @roots.map { |z| laguerre(z) }
> end
> @roots
> end
> end
>
> def irr(a, precision=4)
> p = Polynomial.new(a.reverse)
> p.roots.reject do |r|
> r.image.abs > p.epsilon or r.real <= 0
> end.map do |r|
> sprintf("%.#{precision.to_i}f", (r.real - 1)).to_f
> end
> end
>