discount_to_spot_rates = function(dk){
    ks = 1:length(dk)
    sk = dk^(-1/ks) - 1
    return(sk)
}

## Duplicate Table 4.4 in Example 4.8:
##
year = 1:12
discount = c(0.929, 0.853, 0.776, 0.7, .628, .56, .496, .439, .386, .339, .297, .26)
## plot(discount, type='l')
## grid()

spot = discount_to_spot_rates(discount)

## Consider each bond:
##
bond_1_payments = c(rep(6, 11), 100+6)
PV_1s = discount * bond_1_payments
negPV_1 = year * bond_1_payments * ( 1 + spot )^(-(year+1))

bond_2_payments = c(rep(10, 4), 100+10, rep(0, length(year)-5))
PV_2s = discount * bond_2_payments
negPV_2 = year * bond_2_payments * ( 1 + spot )^(-(year+1) )

## Combine everything to check my numbers:
##
DF = data.frame(year=year, spot=round(spot*100, 2), d=discount,
                B1=bond_1_payments, PV1=PV_1s, negPV1=round(negPV_1, 2),
                B2=bond_2_payments, PV2=PV_2s, negPV2=round(negPV_2, 2) )
print(DF)

PV1 = sum(DF$PV1)
D1 = sum(DF$negPV1)/sum(DF$PV1)
print(sprintf('P_1= %8.3f; negP_1= %8.3f; negDQ_1= %8.3f', PV1, sum(DF$negPV1), D1))

PV2 = sum(DF$PV2)
D2 = sum(DF$negPV2)/sum(DF$PV2)
print(sprintf('P_2= %8.3f; negP_2= %8.3f; negDQ_2= %8.3f', PV2, sum(DF$negPV2), D2))

## The present va1ue of the given obligation stream:
##
obligations = c(500, 900, 600, 500, 100, 100, 100, 50)
obligations = c(obligations, rep(0, length(discount)-length(obligations)))
PV = sum( obligations * discount )
negPV = year * obligations * ( 1 + spot )^(-(year+1))
D = sum(negPV)/PV
print(sprintf('PV= %8.3f; D= %8.3f', PV, D))

A = matrix(data=c(PV1, PV2, PV1*D1/PV, PV2*D2/PV), nrow=2, ncol=2, byrow=T)
#print(A)
b = c(PV, D)
#print(b)
x = solve(A, b)
print(x)