Re: [lmi] Numerics

[Greg Chicares]

On 2016-03-28 22:08, Vadim Zeitlin wrote:
> On Mon, 28 Mar 2016 21:49:38 +0000 Greg Chicares <address@hidden> wrote:
> 
> GC> (Here's the motivation. If a textcontrol is to accept an interest rate in
> GC> a range like [0.03, 0.07], and treat input outside that range as invalid,
> GC> then it would be naive to compare the input directly to either bound: if
> GC> "0.07" is correctly entered and the machine translates that text string to
> GC> 0.07000000000000001 (e.g.), and compares it to a bound represented as, 
> say,
> GC> 0.06999999999999999, then we shouldn't reject it, because it's the user's
> GC> best attempt at specifying the upper-limit interest rate. Therefore, the
> GC> bounds are "adjusted" outward by one ulp.)
> 
>  Unless I'm missing something, the naïve solution should work as long as
> the values are rounded up to the appropriate number of digits. Wouldn't
> this be simpler?

Which values would you round up--input values, or bounds?

If input: there is no explicit maximum number of digits. Rounding may (or may 
not)
be applied to intermediate values calculated from input, yet not to input 
itself.
Sometimes we really do want to enter values to high precision--e.g., if interest
is given as the rounded monthly rate 0.00246627, and we want to enter not its
notional annual equivalent 0.03 but its actual equivalent 0.0300000028... .

If bounds: that's what the present code is intended to do. If (n+1) percent is
represented as 0.0n9999... because that's the closest representable double, then
we want nextafter(0.0n9999...). Calling round(0.0n9999..., 2) is likely to 
return
its argument: not helpful for (n+1) percent when (n+1) is already an integer.

It would be nice if reading "0.07" from a file into the upper bound produced
exactly the same value as reading "0.07" from a text control into an interest
rate, but we have observed that it may not. I can't explain why that happened, 
but
I did need to fix it. With this technique, that problem has not been observed.

> GC> More idealistically, we should probably use integral cents as our basic
> GC> currency unit instead of floating-point dollars rounded to the closest
> GC> approximation to integral hundredths, because in the real world we can
> GC> have exactly seven cents, and (double)(0.07) is not exactly the same.
> 
>  Yes, currency amounts are a classic example of things not to use floating
> point numbers for. I guess there must have been some reasons
> (compatibility?) to do it at some time,

Lack of insight.

> but it would be clearly better to
> avoid it if possible. Please let me know if I should make an issue for this
> too.

If you can think of a tidy strategy for resolving it.

First, let me state the benefits I would hope to gain. Examining one of our
regression-test files, I see values like these:
  BaseDeathBft
  14857345.859999999404
  15767322.6699999999255
  16718005.9499999992549
  17708938.820000000298
They've all been rounded to cents, as the problem domain requires: an insurance
company cannot pay a benefit in fractions of a cent. But 0.01 isn't exactly
representable, so neither are these numbers. Slight changes to the code like
rearranging an expression, or updating the compiler, can cause these calculated
values to change. If some value is calculated as (A-B), rounded down to whole
cents, and originally the values were:
  A = 12.34999
  B =  2.34999
     $10.00 difference
but then B becomes 2.350001, then the difference (A-B) becomes 9.99999, and
rounding it down gives $9.99 instead. Because our calculations are generally
iterative, a one-cent monthly difference becomes twelve cents after a year, and
the change percolates into all later years. That makes acceptance testing hard.
If we could avoid spurious differences like these, we'd save a lot of time.
(And we could make the code simpler by removing kludges like multiplying (A-B)
by (1 + ε) that attempt to "fix" this rounded-subtraction problem.)

In principle, all we have to do is identify which variables, arguments, and
return types are currency, and scale them all by 100, so that currency amounts
could be represented as integral cents rather than fractional dollars. (Scale
them internally only, that is: end users won't enter for $12.34 as "1234", and
customers won't accept "1234¢" on output.)

Glancing at only one header ('account_value.hpp'), I see many things like

  double InforceLivesEoy () const;
  double  SepAcctPaymentAllocation;

where "double" is correct, and many others, e.g.,

  double GetSepAcctAssetsInforce () const;
  void IncrementAVProportionally(double);
  double  SepAcctValueAfterDeduction;

where "double" should really be a currency type.

It wouldn't be too hard to go through all the code and replace "double" with
"currency" as appropriate. But I'm sure I'd make some mistakes. How could I
reliably treat them all, so that I don't break the system?

I was thinking of doing it in steps. First, I could
  typedef double currency;
and s/double/currency/ as seems appropriate. That's safe: even if I misidentify
the type of some quantity, the two types are just synonyms for now. Then...

Here's an idea for checking that the types are all correct: use the compiler as
as a type-enforcement tool. As an intermediate step (not for production release)
we might replace the typedef with a UDT like:

  class currency
  {
    public:
      currency();

      // Probably don't define these, to avoid implicit conversions
      // currency(double);
      // operator=(double);
      // operator double();

      operator+(currency); // ...and other additive operations
      operator*(double);   // ...and other multiplicative operations

    private:
      double value_;
  }

which would have (almost) the same semantics as double, but would be a distinct
type with no implicit conversions. Probably most operations are additive:
  currency = currency + currency - currency;
and most of the rest are multiplicative:
  currency = currency * double / double;

With such a framework, we could even write iostream inserters and extractors
that scale by 100, and run available regression tests to make sure everything
is correct.

Then, for production, we'd push the scaling into the code outside this
"currency" class, and then either go back to
  typedef double currency; // integers stored as double
or use something else like
  typedef int64_t currency;
or even try both and see which runs faster. I think we ultimately want one
of these builtin types for speed and simplicity. We really mustn't make lmi
run any slower. Perhaps we could use expression templates to keep currency
as a class without a speed penalty, but that wouldn't be simple.

Or maybe this idea of temporarily introducing a "currency" class is more
trouble than it's worth?