# 4.4.jl
# 2024-09-24, 動いた
# $Id: 4.4.jl 1.2 2024/12/07 04:17:50 s Exp $

using Dates, Printf

# input
const face_value  = 1000
const cpn_rate = [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]
const ask_price   = [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]
const today       = Date(2011, 11, 05)
const last_cpn = Date(2011, 08, 15)
const next_cpn = Date(2012, 02, 15)

# calculation
const cpn =  face_value .* cpn_rate / 100 / 2
const t = [(Dates.value(Date(2012,2,15)) - Dates.value(today)) / 365, # years
           (Dates.value(Date(2012,8,15)) - Dates.value(today)) / 365,
           (Dates.value(Date(2013,2,15)) - Dates.value(today)) / 365,
           (Dates.value(Date(2013,8,15)) - Dates.value(today)) / 365,
           (Dates.value(Date(2014,2,15)) - Dates.value(today)) / 365,
           (Dates.value(Date(2014,8,15)) - Dates.value(today)) / 365,
           (Dates.value(Date(2015,2,15)) - Dates.value(today)) / 365,
           (Dates.value(Date(2015,8,15)) - Dates.value(today)) / 365,
           (Dates.value(Date(2016,2,15)) - Dates.value(today)) / 365]
const number_of_days_since_last_cpn        =
    Dates.value(today) - Dates.value(last_cpn) # days
const number_of_days_in_current_cpn_period =
    Dates.value(next_cpn) - Dates.value(last_cpn) # days
const accrued_interest = number_of_days_since_last_cpn /
                         number_of_days_in_current_cpn_period .*
                         cpn
const total_price = face_value .* ask_price / 100 .+ accrued_interest

const cf1 = [cpn[ 1]+face_value,0,0,0,0,0,0,0,0]
const cf2 = [cpn[ 2]+face_value,0,0,0,0,0,0,0,0]
const cf3 = [cpn[ 3],cpn[ 3]+face_value,0,0,0,0,0,0,0]
const cf4 = [cpn[ 4],cpn[ 4]+face_value,0,0,0,0,0,0,0]
const cf5 = [cpn[ 5],cpn[ 5],cpn[ 5]+face_value,0,0,0,0,0,0]
const cf6 = [cpn[ 6],cpn[ 6],cpn[ 6]+face_value,0,0,0,0,0,0]
const cf7 = [cpn[ 7],cpn[ 7],cpn[ 7],cpn[ 7]+face_value,
             0,0,0,0,0]
const cf8 = [cpn[ 8],cpn[ 8],cpn[ 8],cpn[ 8]+face_value,
             0,0,0,0,0]
const cf9 = [cpn[ 9],cpn[ 9],cpn[ 9],cpn[ 9],cpn[ 9]+
             face_value,0,0,0,0]
const cf10 = [cpn[10],cpn[10],cpn[10],cpn[10],cpn[10]+
             face_value,0,0,0,0]
const cf11= [cpn[11],cpn[11],cpn[11],cpn[11],cpn[11],
             cpn[11]+face_value,0,0,0]
const cf12= [cpn[12],cpn[12],cpn[12],cpn[12],cpn[12],
             cpn[12]+face_value,0,0,0]
const cf13= [cpn[13],cpn[13],cpn[13],cpn[13],cpn[13],
             cpn[13],cpn[13]+face_value,0,0]
const cf14= [cpn[14],cpn[14],cpn[14],cpn[14],cpn[14],
             cpn[14],cpn[14]+face_value,0,0]
const cf15= [cpn[15],cpn[15],cpn[15],cpn[15],cpn[15],
             cpn[15],cpn[15],cpn[15]+face_value,0]
const cf16= [cpn[16],cpn[16],cpn[16],cpn[16],cpn[16],
             cpn[16],cpn[16],cpn[16]+face_value,0]
const cf17= [cpn[17],cpn[17],cpn[17],cpn[17],cpn[17],
             cpn[17],cpn[17],cpn[17],cpn[17]+face_value]
const cf18= [cpn[18],cpn[18],cpn[18],cpn[18],cpn[18],
             cpn[18],cpn[18],cpn[18],cpn[18]+face_value]

w1 = [cf1;  cf2;  cf3;  cf4;  cf5;  cf6;  cf7;  cf8;  cf9; 
     cf10; cf11; cf12; cf13; cf14; cf15; cf16; cf17; cf18]
w1 = reshape(w1, (9, 18))
w1 = transpose(w1)
const y1 =  total_price

# least suquare method
#        [w1]{beta1} =  {y1}
# [w1]^T [w1]{beta1} =  [w1]^T {y1}
#            {beta1} = ([w1]^T [w1])^-1 [w1]^T]{y1}
A1 = transpose(w1) * w1
b1 = transpose(w1) * y1
beta1 = A1 \ b1
const d = beta1

const r = -log.(d) ./ t
w2 = [map(x -> t[1] ^ x, [0, 1, 2, 3, 4]);
      map(x -> t[2] ^ x, [0, 1, 2, 3, 4]);
      map(x -> t[3] ^ x, [0, 1, 2, 3, 4]);
      map(x -> t[4] ^ x, [0, 1, 2, 3, 4]);
      map(x -> t[5] ^ x, [0, 1, 2, 3, 4]);
      map(x -> t[6] ^ x, [0, 1, 2, 3, 4]);
      map(x -> t[7] ^ x, [0, 1, 2, 3, 4]);
      map(x -> t[8] ^ x, [0, 1, 2, 3, 4]);
      map(x -> t[9] ^ x, [0, 1, 2, 3, 4])]
w2 = reshape(w2, (5, 9))
w2 = transpose(w2)
const y2 = r
const A2 = transpose(w2) * w2
const b2 = transpose(w2) * y2
const beta2 = A2 \ b2

# output
println("4.4.jl")
println("face value  ",  face_value)
println("last_cpn ", last_cpn)
println("today       ",       today)
println("next cpn ", next_cpn)
println("number_of_days_since_last_cpn= ",
         number_of_days_since_last_cpn)
println("number_of_days_in_current_cpn_period= ",
         number_of_days_in_current_cpn_period)
println()
@printf("%5s %10s %10s %10s %15s\n",
    "i", "a_price", "cpn", "a.interest", "total_price")
println(repeat("-", 54))
for i = 1:length(ask_price)
    @printf("%5d %10d %10d %10d %15d\n",
     i, ask_price[i], cpn[i], accrued_interest[i], total_price[i])
end
println(repeat("-", 54)) 
println("a_price = ask price")
println("a.interest = accrued interest")
println()

# print t
@printf("%10s%10s\n", "i", "t(year)")
println(repeat("-", 20)) 
for i = 1:length(t)
    @printf("%10d%10.3f\n", i, t[i])
end
println(repeat("-", 20))
println()

# print y1
@printf("%10s%10s\n", "i", "y1")
println(repeat("-", 20)) 
for i = 1:length(y1)
    @printf("%10d%10d\n", i, y1[i])
end
println(repeat("-", 20))
println()

# print w1
println("w1, matrix of cf(in integer) in row")
@printf("%5s%5s%5s%5s%5s%5s%5s%5s%5s%5s\n",
    'i', "t0","t1", "t2","t3", "t4","t5", "t6", "t7", "t8")
println(repeat("-", 50))
for i = 1:size(w1, 1)
    @printf("%5d", i)
    for j = 1:size(w1, 2)
        @printf("%5d", w1[i, j])
    end
    println()
end
println(repeat("-", 50))
println()

# print beta1
println("estimated beta1, d(t[i])")
@printf("%10s%10s\n", "i", "beta1")
println(repeat("-", 20))
for i = 1:length(beta1)
    @printf("%10d%10.3f\n", i, beta1[i])
end
println(repeat("-", 20))
println()

# print r
println("estimated r(t[i])")
@printf("%10s%10s\n", "i", "r(t[i])")
println(repeat("-", 20))
for i = 1:length(r)
    @printf("%10d%10.4f\n", i, r[i])
end
println(repeat("-", 20))
println()

# print w2
println("w2, coefficent matrix")
@printf("%8s%8s%8s%8s%8s%8s\n", "i", "1", "t", "t^2", "t^3", "t^4")
println(repeat("-", 48))
for i = 1:size(w2, 1)
    @printf("%8d", i)
    for j = 1:size(w2, 2)
        @printf("%8.2f", w2[i, j])
    end
    println()
end
println(repeat("-", 48))
println()

# print beta2
println("beta2, edtimated a[i] for r[t]")
@printf("%10s%15s\n", "i", "a[i]")
println(repeat("-", 25))
for i = 1:length(beta2)
    @printf("%10d%15.5f\n", i, beta2[i])
end
println(repeat("-", 25))
println()

# print r equation estimated
a = beta2
println("estimated r(t)=")
@printf("%.5f + %.5f t + %.5f t^2 + %.5f t^3 + %.5f t^4\n",
    a[1], a[2], a[3], a[4], a[5])
# end