Calculating interest w/ linear algebra

poc.clj

(ns poc
  (:gen-class)
  (:require [clojure.pprint]
            [clojure.string :as str])
  (:import (java.math RoundingMode MathContext)
           (clojure.lang Ratio PersistentVector IPersistentCollection)))

(def ^:dynamic *precision* 10)
(def ^:dynamic *scale* 2)
(def ^:dynamic *rounding* RoundingMode/HALF_EVEN)
(def ^:dynamic *context* (MathContext. *precision* *rounding*))

;;=== scalar math

(defn exp [^BigDecimal b ^Long n]
  (.pow b n *context*))

(defn div [^BigDecimal a ^BigDecimal b]
  (.divide ^BigDecimal a ^BigDecimal b ^MathContext *context*))

(assert (= 3.141443493M (div 10M 3.18325M)))

(defn compounding-factor
  [rate periods]
  {:pre [(not (neg? periods))]}
  (cond (decimal? periods)
        (compounding-factor rate (rationalize periods))

        (ratio? periods)
        (div (dec (exp (inc rate) (.numerator ^Ratio periods)))
             (bigdec (.denominator ^Ratio periods)))

        (integer? periods)
        (dec (exp (inc rate) periods))))

(assert (= 0.03541666667M (compounding-factor 0.425M 1/12)))
(assert (= 0.425760887M (compounding-factor 0.03M 12)))

;;=== linear algebra

(defn transpose' [r m]
  (into (empty m)
        (apply mapcat vector (partition r m))))

(assert (= [1 4
            2 5
            3 6]
           (transpose' 3
                       [1 2 3
                        4 5 6])))

(defn dotp [a b]
  (reduce + (map * a b)))

(assert (= 3 (dotp [1 3 -5] [4 -2 -1])))

(defn dotr [m n a b]
  (into (empty a)
        (map dotp
             (apply interleave (repeat n (partition m a)))
             (apply concat (repeat n (partition m b))))))

(assert (= [3]
           (dotr 3 1
                 [1 3 -5]
                 [4 -2 -1])))

(assert (= [58 64
            139 154]
           (dotr 3 2
                 [1 2 3
                  4 5 6]
                 (transpose' 2 [7 8
                                9 10
                                11 12]))))

(defprotocol LinearAlgebra
  (rank [self])
  (transpose [self])
  (dot [a b]))

;; A generalized matrix type
(deftype Matrix [r v]
  IPersistentCollection
  (seq [self]
    (seq v))
  (count [self]
    (count v))
  (empty [self]
    (Matrix. (second (rank self)) (empty v)))
  (cons [self o]
    (Matrix. r (conj v o)))
  (equiv [self o]
    (.equiv ^IPersistentCollection o v))
  LinearAlgebra
  (rank [self]
    [r (/ (count v) r)])
  (transpose [self]
    (Matrix. (second (rank self))
             (transpose' (first (rank self)) v)))
  (dot [a b]
    (assert (= (first (rank a)) (second (rank b))))
    (dotr (first (rank a)) (second (rank a)) a b)))

(defn M
  [rank array]
  {:pre [(zero? (mod (count array) rank))]}
  (Matrix. rank (vec array)))

(defmethod print-method Matrix [^Matrix o ^java.io.Writer w]
  (->> (.v o)
       (take (- (* 10 (.r o)) 1))
       (partition (.r o) (.r o) (repeat "…"))
       (map (partial str/join \tab))
       (cons (.toString o))
       (str/join \newline)
       (.write w)))

(M 3 (cons 1 (cons 1 (cons 1 (M 3 (range 9))))))

(assert (= (dot (M 3 [1 2 3
                      4 5 6])
                (M 2 [7 9 11 8 10 12]))
           (M 2 [58 64
                 139 154])))

;; Extend vector to always be a row-vector
(extend-type PersistentVector
  LinearAlgebra
  (rank [self]
    [(count self) 1])
  (transpose [self]
    (Matrix. (second (rank self)) self))
  (dot [a b]
    (assert (and (= (first (rank a)) (second (rank b)))
                 (= (second (rank a)) (first (rank b)))))
    (dotr (first (rank a)) (first (rank b)) a b)))

(assert (= [3000]
           (dot [2 1 0]
                (transpose [1000 1000 1000]))))

;;=== financial math

(defn pmt [pv i n]
  (div (* pv i)
       (- 1 (exp (inc i) (- n)))))

(comment
  (pmt 1000M 0.03M 4))

#_(let [c       1000M
        i       0.03M
        n       4
        is      (map (partial compounding-factor i) (range (inc n)))
        iis     (scalar inc is)
        p       (pmt (- c) i n)
        pmts    [0 0 0 0
                 p 0 0 0
                 p p 0 0
                 p p p 0
                 p p p p]
        balance (sum (dot 1 iis [c])
                     (dot 4 iis pmts))
        amort   (scalar (partial + (- c))
                        balance)]
    [balance
     amort])