;;; 4.4.scm
;;; 2016-04-03
;;; $Id: 4.4.scm 1.2 2024/11/02 00:56:51 s Exp $

(use gauche.array)
(use gauche.sequence) ; for 'map-with-index' procedure

(define (main args)
    (define gauss ; 戸川隼人『数値計算』岩波書店, 1991, p.80, ガウスの消去法.
        (lambda (aa bb n )
            (define n1 (+ n 1))
            (define a (make-array (shape 1 n1 1 n1) 0))
            (define b (make-array (shape 1 n1 1  2) 0))
            (define q 0.0)
            (define s 0.0)
            ; 添え字の始まりを1に変更
            (do ((i 1 (+ i 1))) ((> i n) i)
                (do ((j 1 (+ j 1))) ((> j n) j)
                    (array-set! a i j (array-ref aa (- i 1) (- j 1)))
                )
            )
            (do ((i 1 (+ i 1))) ((> i n) i)
                (array-set! b i 1 (array-ref bb (- i 1) 0))
            )
            (set! q 0.0)
            (set! s 0.0)
            ; 前進消去
            (do ((k 1 (+ k 1))) ((> k (- n 1)) k) ; do k=1,n-1
                (do ((i (+ k 1) (+ i 1))) ((> i n) i)
                    (set! q (/ (array-ref a i k) (array-ref a k k)))
                    (do ((j k (+ j 1))) ((> j n) j) ; j=k,n
                        (array-set! a i j (- (array-ref a i j) (* q (array-ref a k j))))
                    )
                    (array-set! b i 1 (- (array-ref b i 1) (* q (array-ref b k 1))))
                )
            )
            ; 後退代入
            (array-set! x1 n 1(/ (array-ref b n 1) (array-ref a n n)))
            (do ((k (- n 1) (- k 1))) ((< k 1) k)  ; k=n-1,1,-1
                (set! s (array-ref b k 1))
                (do ((j (+ k 1) (+ j 1))) ((> j n) j) ; j=k+1,n
                    (set! s (- s (* (array-ref a k j) (array-ref x1 j 1))))
                )
                (array-set! x1 k 1 (/ s (array-ref a k k)))
            )
        )
    )
    (define gauss2 ; 戸川隼人『数値計算』岩波書店, 1991, p.80, ガウスの消去法.
        (lambda (aa bb n )
            (define n1 (+ n 1))
            (define a (make-array (shape 1 n1 1 n1) 0))
            (define b (make-array (shape 1 n1 1  2) 0))
            (define q 0.0)
            (define s 0.0)
            ; 添え字の始まりを1に変更
            (do ((i 1 (+ i 1))) ((> i n) i)
                (do ((j 1 (+ j 1))) ((> j n) j)
                    (array-set! a i j (array-ref aa (- i 1) (- j 1)))
                )
            )
            (do ((i 1 (+ i 1))) ((> i n) i)
                (array-set! b i 1 (array-ref bb (- i 1) 0))
            )
            ; 前進消去
            (do ((k 1 (+ k 1))) ((> k (- n 1)) k) ; do k=1,n-1
                (do ((i (+ k 1) (+ i 1))) ((> i n) i)
                    (set! q (/ (array-ref a i k) (array-ref a k k)))
                    (do ((j k (+ j 1))) ((> j n) j) ; j=k,n
                        (array-set! a i j (- (array-ref a i j) (* q (array-ref a k j))))
                    )
                    (array-set! b i 1 (- (array-ref b i 1) (* q (array-ref b k 1))))
                )
            )
            ; 後退代入
            (array-set! x2 n 1 (/ (array-ref b n 1) (array-ref a n n)))
            (do ((k (- n 1) (- k 1))) ((< k 1) k) ; k=n-1,1,-1
                (set! s (array-ref b k 1))
                (do ((j (+ k 1) (+ j 1))) ((> j n) j) ; j=k+1,n
                    (set! s (- s (* (array-ref a k j) (array-ref x2 j 1))))
                )
                (array-set! x2 k 1 (/ s (array-ref a k k)))
            )
        )
    )

    (define mjd (lambda (y m d) ; 修正ユリウス日, 参考文献 wikipedia
        (define mjd2 (lambda (yy mm dd) 
                             (+ (floor (* 365.25 yy)) 
                                (- (floor (/ yy 400.0)) (floor (/ yy 100.0)))
                                (floor (* 30.59 (- mm 2)))
                                (- d 678912))))
        (cond ([= m 1] (let ((yy (- y 1)) (mm 13) (dd d)) (mjd2 yy mm dd)))
              ([= m 2] (let ((yy (- y 1)) (mm 14) (dd d)) (mjd2 yy mm dd)))
              (else    (let ((yy    y   ) (mm  m) (dd d)) (mjd2 yy mm dd)))
    )))
    
    (define (fix n x) (/ (round (* x (expt 10 n))) (expt 10 n)))

    ;; input
    (define face_value 1000)
    (define coupon_rate 
        (vector (+  6 (/ 5 8)) (+  9 (/ 1 8)) (+ 7 (/ 7 8))
                (+  8 (/ 1 4)) (+  8 (/ 1 4)) (+ 8 (/ 3 8))
                    8          (+  8 (/ 3 4)) (+ 6 (/ 7 8))
                (+  8 (/ 7 8)) (+  6 (/ 7 8)) (+ 8 (/ 5 8))
                (+  7 (/ 3 4)) (+ 11 (/ 1 4)) (+ 8 (/ 1 2))
                (+ 10 (/ 1 2)) (+  7 (/ 7 8)) (+ 8 (/ 7 8))))
    (define ask_price 
        (vector (+ 100 (/  0 32)) (+ 100 (/ 22 32)) (+ 100 (/ 24 32))
              (+ 101 (/  1 32)) (+ 101 (/  7 32)) (+ 101 (/ 12 32))
              (+ 100 (/ 26 32)) (+ 102 (/  1 32)) (+  98 (/  5 32))
              (+ 102 (/  9 32)) (+  97 (/ 13 32)) (+ 101 (/ 23 32))
              (+  99 (/  5 32)) (+ 109 (/  4 32)) (+ 101 (/ 13 32))
              (+ 107 (/ 27 32)) (+  99 (/ 13 32)) (+ 103 (/  0 32))))
    (define today (mjd 2011 11 5))
;(define today (date->modified-julian-day (make-date 0 0 0 0 5 11 2011 (date-zone-offset (current-date)))))
    (define last_coupon (mjd 2011 8 15))
    (define next_coupon (mjd 2012 2 15))
    (define coupon (map (lambda (x)(/ (* face_value x) 100 2)) coupon_rate))

    (define t2 (vector (- (mjd 2012 2 15) today) (- (mjd 2012 8 15) today) 
                     (- (mjd 2013 2 15) today) (- (mjd 2013 8 15) today)
                     (- (mjd 2014 2 15) today) (- (mjd 2014 8 15) today)
                     (- (mjd 2015 2 15) today) (- (mjd 2015 8 15) today)
                     (- (mjd 2016 2 15) today)))
    (define t (map (lambda (x) (/ x 365)) t2))

    (define number_of_days_since_last_coupon        (- today        last_coupon))
    (define number_of_days_in_current_coupon_period (- next_coupon  last_coupon))
    (define accrued_interest 
        (map (lambda (x) (* (/ number_of_days_since_last_coupon 
                               number_of_days_in_current_coupon_period)
                             x)) coupon))
    (define total_price 
        (map (lambda (x y) (+ (/ (* face_value  x) 100) y))
                               ask_price accrued_interest))

    (define cf0 (vector (+ (ref coupon 0) face_value) 0 0 0 0 0 0 0 0))
    (define cf1 (vector (+ (ref coupon 1) face_value) 0 0 0 0 0 0 0 0))
    (define cf2 (vector (ref coupon 2) (+ (ref coupon 2) face_value) 0 0 0 0 0 0 0))
    (define cf3 (vector (ref coupon 3) (+ (ref coupon 3) face_value) 0 0 0 0 0 0 0))
    (define cf4 (vector (ref coupon 4) (ref coupon 4) (+ (ref coupon 4) face_value)
                    0 0 0 0 0 0))
    (define cf5 (vector (ref coupon 5) (ref coupon 5) (+ (ref coupon 5) face_value)
                    0 0 0 0 0 0))
    (define cf6 (vector (ref coupon 6) (ref coupon 6) (ref coupon 6)
                  (+ (ref coupon 6) face_value) 0 0 0 0 0))
    (define cf7 (vector (ref coupon 7) (ref coupon 7) (ref coupon 7)
                  (+ (ref coupon 7) face_value) 0 0 0 0 0))
    (define cf8 (vector [ref coupon 8] (ref coupon 8) (ref coupon 8) (ref coupon 8)
                  (+ (ref coupon 8) face_value) 0 0 0 0))
    (define cf9 (vector [ref coupon 9] [ref coupon 9] [ref coupon 9] [ref coupon 9]
                  (+ (ref coupon 9) face_value) 0 0 0 0))
    (define cf10 (vector [ref coupon 10] [ref coupon 10] [ref coupon 10]
                   [ref coupon 10] [ref coupon 10] (+ (ref coupon 10) face_value)
                        0 0 0))
    (define cf11 (vector [ref coupon 11] [ref coupon 11] [ref coupon 11]
                   [ref coupon 11] [ref coupon 11] (+ (ref coupon 11) face_value)
                    0 0 0))
    (define cf12 (vector [ref coupon 12] [ref coupon 12] [ref coupon 12]
                   [ref coupon 12] [ref coupon 12] [ref coupon 12]
                   (+ (ref coupon 12) face_value) 0 0))
    (define cf13 (vector [ref coupon 13] [ref coupon 13] [ref coupon 13]
                   [ref coupon 13] [ref coupon 13] [ref coupon 13]
                   (+ [ref coupon 13] face_value) 0 0))
    (define cf14 (vector [ref coupon 14] [ref coupon 14] [ref coupon 14]
                   [ref coupon 14] [ref coupon 14] [ref coupon 14]
                   [ref coupon 14] (+ (ref coupon 14) face_value) 0))
    (define cf15 (vector [ref coupon 15] [ref coupon 15] [ref coupon 15]
                   [ref coupon 15] [ref coupon 15] [ref coupon 15]
                   [ref coupon 15] [+ [ref coupon 15] face_value] 0))
    (define cf16 (vector [ref coupon 16] [ref coupon 16] [ref coupon 16]
                   [ref coupon 16] [ref coupon 16] [ref coupon 16]
                   [ref coupon 16] [ref coupon 16] (+ [ref coupon 16] face_value)))
    (define cf17 (vector [ref coupon 17] [ref coupon 17] [ref coupon 17]
                   [ref coupon 17] [ref coupon 17] [ref coupon 17]
                   [ref coupon 17] [ref coupon 17] [+ [ref coupon 17] face_value]))

    (define y1 (make-f64array (shape 1 19 1  2) 0.0))
    (define y2 (make-f64array (shape 0  9 0  1) 0.0))
    (define w1 (make-f64array (shape 1 19 1 10) 0.0))
    (define w2 (make-f64array (shape 0  9 0  5) 0.0))
    (define x1 (make-f64array (shape 1 10 1  2) 0.0))
    (define x2 (make-f64array (shape 1  6 1  2) 0.0))
    (define d  (list->vector (iota 9)))
    (define r  (list->vector (iota 9)))
    (define a  (vector 1. 2. 3. 4. 5.))

    ; y1
    (do ((i 1 (+ i 1)))
        ((> i 18) i)
        (array-set! y1  i 1  (ref  total_price (- i 1)))
    )
    ; w1
    (do ((j 1 (+ j 1))) (( > j 9) j)
        (array-set! w1  1 j  (ref  cf0 (- j 1)))
        (array-set! w1  2 j  (ref  cf1 (- j 1)))
        (array-set! w1  3 j  (ref  cf2 (- j 1))) 
        (array-set! w1  4 j  (ref  cf3 (- j 1)))
        (array-set! w1  5 j  (ref  cf4 (- j 1))) 
        (array-set! w1  6 j  (ref  cf5 (- j 1)))
        (array-set! w1  7 j  (ref  cf6 (- j 1))) 
        (array-set! w1  8 j  (ref  cf7 (- j 1)))
        (array-set! w1  9 j  (ref  cf8 (- j 1))) 
        (array-set! w1 10 j  (ref  cf9 (- j 1)))
        (array-set! w1 11 j  (ref cf10 (- j 1))) 
        (array-set! w1 12 j  (ref cf11 (- j 1)))
        (array-set! w1 13 j  (ref cf12 (- j 1))) 
        (array-set! w1 14 j  (ref cf13 (- j 1)))
        (array-set! w1 15 j  (ref cf14 (- j 1))) 
        (array-set! w1 16 j  (ref cf15 (- j 1)))
        (array-set! w1 17 j  (ref cf16 (- j 1))) 
        (array-set! w1 18 j  (ref cf17 (- j 1)))
    )
    ; beta1
    (gauss (array-mul (array-transpose w1) w1) 
           (array-mul (array-transpose w1) y1) 9)
    (set! d (array->vector x1))
    (set! r (map (lambda (x y) (/ (- (log x)) y)) d t))
    ; w2
    (do ((i 0 (+ i 1))) ((> i 8) i)
        (do ((j 0 (+ j 1))) ((> j 4) j)
            (array-set! w2 i j 
                (ref (map (lambda (x) (expt (ref t i) x))(iota 5)) j))
        )
    )
    ; y2
    (do ((i 0 (+ i 1))) 
        (( > i 8) i)
        (array-set! y2  i 0  (ref r i))
    )
    ; beta2
    (gauss2 (array-mul (array-transpose w2) w2) 
            (array-mul (array-transpose w2) y2) 5)
    (set! a (array->vector x2))

    ;; output
    (print "4.4.scm")
    (format #t "~10@a~10@a~10@a~10@a~10@a\n" 'i 'a_price 'coupon 'a.int 'total_pri)
    (print "--------------------------------------------------")
    (for-each-with-index (^(i x1 x2 x3 x4) 
        (format #t "~10d~10,1f~10,2f~10,2f~10,2f\n" i (exact->inexact x1) (exact->inexact x2) x3 x4))
         ask_price coupon accrued_interest total_price)
    (print "--------------------------------------------------")
    (print "a_price = ask price")
    (print "a.int = accrued interest")(newline)

    (format #t "~10@a~10@a\n" 'i 't)
    (print "--------------------")
    (for-each-with-index 
        (lambda (i x) (format #t "~10d~10,3f\n" i x)) t)
    (print "--------------------")(newline)

    (format #t "~10@a~10@a\n" 'i "y1")
    (print "--------------------")
    (for-each-with-index 
        (lambda (i x) (format #t "~10d~10d\n" i (x->integer x))) (array->vector y1))
    (print "--------------------")(newline)

    (print "w1, a matrix of cf(in integer)")
    (print "---------------------------------------------")
    (do ((i 1 (+ i 1))) ((>= i 19) i)
        (do ((j 1 (+ j 1))) ((>= j 10) j)
            (format #t "~5d" (x->integer (array-ref w1 i j)))
        )
        (newline)
    )
    (print "---------------------------------------------")(newline)

    (print "an estimated beta1 (=d(t[i]))")
    (format #t "~10@a~10@a\n" 'i "beta[i]")
    (print "--------------------")
    (for-each-with-index (lambda (i x) (format #t "~10d~10,3f\n" i x))
         (array->vector x1))
    (print "--------------------")(newline)

    (print "an estimated r")
    (format #t "~10@a~10@a\n" 'i "r[i]")
    (print "--------------------")
    (for-each-with-index (lambda (i x) (format #t "~10d~10,4f\n" i x)) r)
    (print "--------------------")(newline)

    (print "w2, a coefficient matrix")
    (format #t "~10@a~10@a~10@a~10@a~10@a\n" '1 't 't^2 't^3 't^4)
    (print "--------------------------------------------------")
    (do ((i 0 (+ i 1))) ((> i 8) i)
        (do ((j 0 (+ j 1))) ((> j 4) j)
            (format #t "~10,3f" (array-ref w2 i j))
        )
        (newline)
    )
    (print "--------------------------------------------------")
    (newline)
    
    (print "beta2, an estimated a[i] for r(t)")
    (format #t "~10@a~20@a\n" 'i 'beta2)
    (print "------------------------------")
    (for-each-with-index  (lambda (i x) (format #t "~10d~20,5f\n" i x))
        (array->vector x2)
    )
    (print "------------------------------")(newline)

    (print "r(t)=")
    (format #t "~7,5f + ~7,5f t + ~7,5f t^2 + ~8,5f t^3 + ~7,5f t^4" 
        (vector-ref a 0)(vector-ref a 1)(vector-ref a 2)(vector-ref a 3)
        (vector-ref a 4)
    )
0)
; end