Let x1,x2,x3,x4,x5,x6,x7,x8,x9,x10 represent the initial investment in each of the ten bonds. The objective is to minimize the initial cost and matching the liability cash flows for the coming 8 years.
\(min \ 102x_1 + 99x_2 + 101x_3 + 98x_4 + 98x_5 + 104x_6 + 100x_7 + 101x_8 + 102x_9 + 94x_{10}\)
\(s.t.\)
\(Year\ 1:\ (100+5)x_1 + 3.5x_2 + 5x_3 + 3.5x_4 + 4x_5 + 9x_6 + 6x_7 + 8x_8 + 9x_9 + 7x_{10}\ =\ 12,000\)
\(Year\ 2:\ (100+3.5)x_2 + (100+5)x_3 + 3.5x_4 + 4x_5 + 9x_6 + 6x_7 + 8x_8 + 9x_9 + 7x_{10}\ =\ 18,000\)
\(Year\ 3:\ (100+3.5)x_4 + 4x_5 + 9x_6 + 6x_7 + 8x_8 + 9x_9 + 7x_{10}\ =\ 20,000\)
\(Year\ 4:\ (100+4)x_5 + 9x_6 + 6x_7 + 8x_8 + 9x_9 + 7x_{10}\ =\ 20,000\)
\(Year\ 5:\ (100+9)x_6 + (100+6)x_7 + 8x_8 + 9x_9 + 7x_{10}\ =\ 16,000\)
\(Year\ 6:\ (100+8)x_8 + 9x_9 + 7x_{10}\ =\ 15,000\)
\(Year\ 7:\ (100+9)x_9 + 7x_{10}\ =\ 12,000\)
\(Year\ 8:\ (100+7)x_{10}\ =\ 10,000\)
library(lpSolveAPI)
my.lp <- make.lp(8,10) # 8 Constraints, 10 decision variables
print (my.lp)
## Model name:
## a linear program with 10 decision variables and 8 constraints
Colnames=c('x1','x2','x3','x4','x5','x6','x7','x8','x9','x10')
Rownames=c('Y1','Y2','Y3','Y4','Y5','Y6','Y7','Y8')
dimnames(my.lp)<-list(Rownames,Colnames)
add.constraint(my.lp,c(105,3.5,5,3.5,4,9,6,8,9,7), "=",12000) #Y1
add.constraint(my.lp,c(0,103.5,105,3.5,4,9,6,8,9,7), "=",18000) #Y2
add.constraint(my.lp,c(0,0,0,103.5,4,9,6,8,9,7), "=",20000) #Y3
add.constraint(my.lp,c(0,0,0,0,104,9,6,8,9,7), "=",20000) #Y4
add.constraint(my.lp,c(0,0,0,0,0,109,106,8,9,7), "=",16000) #Y5
add.constraint(my.lp,c(0,0,0,0,0,0,0,108,9,7), "=",15000) #Y6
add.constraint(my.lp,c(0,0,0,0,0,0,0,0,109,7), "=",12000) #Y7
add.constraint(my.lp,c(0,0,0,0,0,0,0,0,0,107), "=",10000) #Y8
set.objfn(my.lp, c(102,99,101,98,98,104,100,101,102,94))
lp.control(my.lp,sense='min')
## $anti.degen
## [1] "none"
##
## $basis.crash
## [1] "none"
##
## $bb.depthlimit
## [1] -50
##
## $bb.floorfirst
## [1] "automatic"
##
## $bb.rule
## [1] "pseudononint" "greedy" "dynamic" "rcostfixing"
##
## $break.at.first
## [1] FALSE
##
## $break.at.value
## [1] -1e+30
##
## $epsilon
## epsb epsd epsel epsint epsperturb epspivot
## 1e-10 1e-09 1e-12 1e-07 1e-05 2e-07
##
## $improve
## [1] "dualfeas" "thetagap"
##
## $infinite
## [1] 1e+30
##
## $maxpivot
## [1] 250
##
## $mip.gap
## absolute relative
## 1e-11 1e-11
##
## $negrange
## [1] -1e+06
##
## $obj.in.basis
## [1] TRUE
##
## $pivoting
## [1] "devex" "adaptive"
##
## $presolve
## [1] "none"
##
## $scalelimit
## [1] 5
##
## $scaling
## [1] "geometric" "equilibrate" "integers"
##
## $sense
## [1] "minimize"
##
## $simplextype
## [1] "dual" "primal"
##
## $timeout
## [1] 0
##
## $verbose
## [1] "neutral"
write.lp(my.lp,'Project1.lp',type='lp')
solve(my.lp)
## [1] 0
print (get.objective(my.lp))
## [1] 93944.5
print (get.variables(my.lp))
## [1] 62.13613 0.00000 125.24293 151.50508 156.80776 123.08007 0.00000
## [8] 124.15727 104.08986 93.45794
Testing of our function dedicate_g3(P,C,M,L):
P = c(102,99,101,98,98,104,100,101,102,94)
C = c(5,3.5,5,3.5,4,9,6,8,9,7)
M = c(1,2,2,3,4,5,5,6,7,8)
L = c(12000,18000,20000,20000,16000,15000,12000,10000)
Portfolio<-dedicate_g3(P,C,M,L)
print (get.objective(Portfolio))
## [1] 93944.5
print (get.variables(Portfolio))
## [1] 62.13613 0.00000 125.24293 151.50508 156.80776 123.08007 0.00000
## [8] 124.15727 104.08986 93.45794
Get US Treasury Bonds quotes from http://www.wsj.com/mdc/public/page/2_3020-treasury.html
* STARTING DATE: pulled on Feb 8, 2016, which should be treated as the STARTING DATE.
* MATURITY DATES: chosen to match the date of liability stream. For some of the maturity dates there are multiple bonds available, we just randomly pick one bond for each maturity date.
* PRICES: calculated the average of bid and ask price (\(price\ =\ (bid\ +\ ask)\ /\ 2\)) as the price of a bond, which is a widely used approach in academics.
raw_quotes = read.csv("US Treasury Bonds.csv")
print (raw_quotes)
## Maturity Coupon Bid Asked Chg Asked.yield
## 1 6/30/2016 0.500 100.0625 100.0781 0.0234 0.299
## 2 6/30/2016 1.500 100.4609 100.4766 0.0391 0.277
## 3 6/30/2016 3.250 101.1328 101.1484 0.0078 0.302
## 4 12/31/2016 0.625 100.0000 100.0391 unch. 0.605
## 5 12/31/2016 0.875 100.2813 100.2969 0.0625 0.540
## 6 12/31/2016 3.250 102.3672 102.3828 0.0469 0.562
## 7 6/30/2017 0.625 100.0078 100.0234 0.0938 0.608
## 8 6/30/2017 0.750 100.1719 100.1875 0.0859 0.614
## 9 6/30/2017 2.500 102.5859 102.6016 0.0938 0.617
## 10 12/31/2017 0.750 100.1328 100.1484 0.1094 0.671
## 11 12/31/2017 1.000 100.6172 100.6328 0.1016 0.663
## 12 12/31/2017 2.750 103.8750 103.8906 0.1406 0.675
## 13 6/30/2018 1.375 101.5156 101.5313 0.1641 0.728
## 14 6/30/2018 2.375 103.9375 103.9531 0.1797 0.704
## 15 12/31/2018 1.375 101.5781 101.5938 0.2188 0.816
## 16 12/31/2018 1.500 101.9453 101.9609 0.2422 0.812
## 17 6/30/2019 1.000 100.2656 100.2813 0.2656 0.916
## 18 6/30/2019 1.625 102.3359 102.3516 0.2734 0.919
## 19 12/31/2019 1.125 100.3750 100.3906 0.3359 1.022
## 20 12/31/2019 1.625 102.2656 102.2813 0.3047 1.025
## 21 6/30/2020 1.625 102.2109 102.2266 0.3828 1.104
## 22 6/30/2020 1.875 103.3125 103.3281 0.4063 1.096
## 23 12/31/2020 1.750 102.8125 102.8281 0.4453 1.153
## 24 12/31/2020 2.375 105.6328 105.6484 0.4609 1.183
## 25 6/30/2021 2.125 104.4688 104.4844 0.4688 1.262
## 26 12/31/2021 2.125 104.4297 104.4453 0.5938 1.338
P4 = (raw_quotes$Bid + raw_quotes$Asked) / 2
C4 = raw_quotes$Coupon
M4 = c(1,1,1,
2,2,2,
3,3,3,
4,4,4,
5,5,
6,6,
7,7,
8,8,
9,9,
10,10,
11,
12)
L4 = c(9,9,10,10,6,6,9,9,10,10,5,3) * 1000000
Portfolio_4<-dedicate_g3(P4,C4,M4,L4)
print (get.objective(Portfolio_4))
## [1] 87582418
print (get.variables(Portfolio_4))
## [1] 0.00 0.00 70495.75 0.00 0.00 72786.86 0.00
## [8] 0.00 85152.43 0.00 0.00 87281.24 0.00 49681.48
## [15] 0.00 50861.41 0.00 81624.33 0.00 82950.73 0.00
## [22] 94298.68 0.00 96066.78 48348.36 29375.76
## $duals
## [1] 0.9795700 0.9606915 0.9535912 0.9335738 0.9265462 0.9342083 0.9161176
## [8] 0.9007770 0.8760568 0.8374643 0.8312077 0.8135292 1.6235174 1.0423975
## [15] 0.0000000 2.7377364 2.5222210 0.0000000 2.8478237 2.6501922 0.0000000
## [22] 3.9126529 3.4401963 0.0000000 2.3321226 0.0000000 0.3438726 0.0000000
## [29] 2.0573866 0.0000000 1.8618878 0.0000000 0.9937331 0.0000000 2.9413229
## [36] 0.0000000 0.0000000 0.0000000
##
## $dualsfrom
## [1] 1.721314e+06 1.484757e+06 1.271876e+06 1.031852e+06 9.138590e+05
## [6] 8.375669e+05 7.049273e+05 5.701324e+05 3.933224e+05 1.651638e+05
## [11] 6.242350e+04 0.000000e+00 -1.000000e+30 -1.000000e+30 -1.000000e+30
## [16] -2.862950e+06 -3.164313e+06 -1.000000e+30 -4.108333e+06 -4.401785e+06
## [21] -1.000000e+30 -3.957480e+06 -4.522834e+06 -1.000000e+30 -8.102940e+06
## [26] -1.000000e+30 -4.129947e+07 -1.000000e+30 -8.394116e+06 -1.000000e+30
## [31] -1.066315e+07 -1.000000e+30 -2.172617e+07 -1.000000e+30 -8.896867e+06
## [36] -1.000000e+30 -1.000000e+30 -1.000000e+30
##
## $dualstill
## [1] 1.000000e+30 2.402382e+08 3.181250e+08 2.978167e+08 3.471764e+08
## [6] 3.501622e+08 3.318506e+08 3.370969e+08 2.996823e+08 2.441281e+08
## [11] 2.722331e+08 2.402941e+08 7.242473e+04 7.171119e+04 1.000000e+30
## [16] 7.468565e+04 7.450055e+04 1.000000e+30 8.673912e+04 8.663150e+04
## [21] 1.000000e+30 8.901387e+04 8.879354e+04 1.000000e+30 5.017155e+04
## [26] 1.000000e+30 5.092412e+04 1.000000e+30 8.212943e+04 1.000000e+30
## [31] 8.336087e+04 1.000000e+30 9.453065e+04 1.000000e+30 9.665687e+04
## [36] 1.000000e+30 1.000000e+30 1.000000e+30
The first 12 duals correspond to the 12 constraints we have, and the rest 26 correspond to the non-negative constraints.
The first 12 duals can be interpreted as follows: if liability in period i is reduced by $1, the intial investment in period 0 will reduce by \(dual amount\). For example, for the first dual, if liability in year 1 reduces by $1, our initial payout would reduce by $0.97.