Maximum Smoothness Forward Rates by Inverse Matrix using R

[This article was first published on K & L Fintech Modeling, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

This post solves Adams and Deventer (1994) maximum smoothness forward rate curve as the matrix inversion of equality constrained quadratic programming problems. We implement R code for this model and elaborate on how to get monthly forward and spot rates from the calibrated parameters.



Maximum Smoothness Forward Rates



Using Adams and Deventer (1994, revised 2010) approach, output of this post is the following figure describing maximum smoothness forward rate curve with spot curve.

Maximum Smoothness Forward Rates by R code

At first, details of maximum smoothness forward rate curve and the solution by using an optimization are explained in the previous post below.


We can solve Adams and Deventer (1994) maximum smoothness forward rate curve by matrix inversion. This solution is possible due to the fact that this problem can be formulated as the equality constrained quadratic programming problem.


Quadratic Programming with Linear Equality Constraints


Equality-constrained quadratic program (QP) is a QP where only linear equality constraints are present with quadratic terms absent.
minxf(x)=12xTQx+cTxs.t.Ax=b

The first-order condition (FOC) for this problem constitutes the following linear system, which is derived from setting up the Lagrangian.
[QATA0]KKT matrix[xλ]=[cb]
Here, matrix in LHS is called a KKT(Karush–Kuhn–Tucker) matrix. We can solve the above linear system by using the inverse matrix of the LHS coefficient matrix. The derivation of it is as follows.
Using Lagrange multiplier (λ), the optimization problems can be converted to the following Lagrangian. f(x,λ)=12xTQx+cTx+λT(Axb)

The FOC with respect to x and λ is
f(x,λ)x=Qx+c+ATλ=0f(x,λ)λ=Axb=0

This two set of equations are reduced to the above vector-matrix representation and we use the inverse matrix of A to get an vector of optimal decision variables x.
X=[xλ]=[QATA0]1[cb]


Objecive Function


As can be seen in previous blog, these quadratic terms can be summarized as the following vector-matrix notation (Δtni=tnitni1).
titi1[f(t)]2dt=xTiQixi
where xi=[aibicidiei],Qi=[1445Δt5i18Δt4i8Δt3i0018Δt4i12Δt3i6Δt2i008Δt3i6Δt2i4Δti000000000000]

Therefore, the objective function is the sum of i-th objective functions.
xTQx
where x=[x1x2xm1xm],Q=[Q10000Q20000Qm10000Qm]

In this way, concatenation of the objective functions in each segment consist of Q matrix in the above QP.


Linear Equality Constraints


As can be seen also in previous blog, for Δki+1=ki+1ki,Δtni=tnitni1, linear equality constraints are as follows.


1) continuity or smoothness at the joints (internal knots)0=Δai+1t4+Δbi+1t3+Δci+1t2+Δdi+1t+Δei+10=4Δai+1t3+3Δbi+1t2+2Δci+1t+dΔi+10=12Δai+1t2+6Δbi+1t+2Δci+10=24Δai+1t+6Δbi+1

2) market price fitting at the observed pointslog[PiPi1]=15aiΔt5i+14biΔt4i+13ciΔt3i+12diΔt2i+eiΔti

2) Additional constraints for smoothness
f(0)=r(0),f(0)=0f(T)=0,f(T)=0

These linear equality constraints consist of A matrix in the above QP.


Fitting Yield Curve as a Quadratic Programming with Linear Equality Constraints


In our problem, the objective function consists of quadratic terms and the linear constraints are equality restrictions. Therefore, this can be represented as a equality-constrained QP.

[QATA0][xλ]=[cb]
Using Adams and Deventer (2010) data, let’s construct left hand side KKT coefficient matrix and right hand side vector. The former consists of Q, A and AT and the latter is made from b and c. Of course, the unknown X vector contains 25 coefficients of forward rate equations.

a1, b1, c1, d1, e1, a2, b2, c2, d2, e2, a3, b3, c3, d3, e3, a4, b4, c4, d4, e4, a5, b5, c5, d5, e5

Let’s construct three sub matrices : Q, A and AT.

Q matrix is of the next form. For example, entry of (1,1) is 0.028125 = (144/5)*(0.25-0)^5. Values of remaining entries are calculated by using the objective function and constraints.
Maximum Smoothness Forward Rates by using R code

A matrix is of the following form.
Maximum Smoothness Forward Rates by using R code

Its minus transpose, AT is
Maximum Smoothness Forward Rates by Inverse Matrix using R

As vector c is a zero vector and vector b has some values which are derived from bond price matching conditions, the RHS vector consists of the following numbers.

0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0.012 , 0.033125 , 0.12 , 0.0975 , 0.3875 , 0.04 , 0 , 0 , 0


Finally, we can construct LHS coefficient matrix and RHS vector as follows.
Maximum Smoothness Forward Rates by Inverse Matrix using R

As we know value of each entry, entries that have numerical value are colored black, blue and red. In LHS matrix, we display empty cells as zero values (0).


R code


The following R code implements the maximum smoothness forward curve by using an inverse matrix.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
#========================================================#
# Quantitative ALM, Financial Econometrics & Derivatives 
# ML/DL using R, Python, Tensorflow by Sang-Heon Lee 
#
# https://kiandlee.blogspot.com
#——————————————————–#
# Adams and Deventer Maximum Smoothness Forward Curve
# using the inverse matrix
#========================================================#
 
graphics.off()  # clear all graphs
rm(list = ls()) # remove all files from your workspace
 
# Input : market zero rate, maturity
df.mkt < data.frame(
    mat = c(0.2513510),
    zrc = c(4.754.55.55.256.5)/100
)
 
# number of maturity
< length(df.mkt$mat)
    
# add 0-maturity zero rate (assumption)
#df <- rbind(c(0, df.mkt$zrc[1]), df.mkt)
df < rbind(c(00.04), df.mkt)
    
# discount factor
df$DF < with(df, exp(zrc*mat))
    
# -ln(P(t(i)/t(i-1)))
df$mln < c(NA,log(df$DF[1:n+1]/df$DF[1:n]))
    
# ti^n
df$t5 < df$mat^5
df$t4 < df$mat^4
df$t3 < df$mat^3
df$t2 < df$mat^2
df$t1 < df$mat^1
df$t0 < 1
    
# dti = ti^n-(ti-1)^n
df$dt5 < c(NA,df$t5[1:n+1 df$t5[1:n])
df$dt4 < c(NA,df$t4[1:n+1 df$t4[1:n])
df$dt3 < c(NA,df$t3[1:n+1 df$t3[1:n])
df$dt2 < c(NA,df$t2[1:n+1 df$t2[1:n])
df$dt1 < c(NA,df$t1[1:n+1 df$t1[1:n])
    
    
# construction linear system
mQ < mA < matrix(0, nrow = n*n, ncol = n*n)
vC < vB < rep(0,n*n)
    
# Objective function
for(r in 1:n) {
    mQ[((r1)*n+1):(r*n2),((r1)*n+1):(r*n2)] < 
        matrix(with(df[r+1,], 
               c(144/5*dt5, 18*dt4, 8*dt3,
                    18*dt4, 12*dt3, 6*dt2,
                     8*dt3,  6*dt2, 4*dt1)),3,3)
}
 
# Smoothness Constraints : f, f’, f”, f”’
= 1for(t in 1:(n1)) {
    mA[(r1)*(n1)+t,(1+(t1)*n):(t*nr+1)] < 
        with(df[t+1,],  c(t4, t3, t2, t1, t0))
    mA[(r1)*(n1)+t,(1+t*n):((t+1)*nr+1)] < 
        with(df[t+1,], c(t4, t3, t2, t1, t0))
}
    
= 2for(t in 1:(n1)) {
    mA[(r1)*(n1)+t,(1+(t1)*n):(t*nr+1)] < 
        with(df[t+1,],  c(4*t3, 3*t2, 2*t1, t0))
    mA[(r1)*(n1)+t,(1+t*n):((t+1)*nr+1)] < 
        with(df[t+1,], c(4*t3, 3*t2, 2*t1, t0))
}
    
= 3for(t in 1:(n1)) {
    mA[(r1)*(n1)+t,(1+(t1)*n):(t*nr+1)] < 
        with(df[t+1,],  c(12*t2, 6*t1, 2*t0))
    mA[(r1)*(n1)+t,(1+t*n):((t+1)*nr+1)] < 
        with(df[t+1,], c(12*t2, 6*t1, 2*t0))
}
    
= 4for(t in 1:(n1)) {
    mA[(r1)*(n1)+t,(1+(t1)*n):(t*nr+1)] < 
        with(df[t+1,],  c(24*t1, 6*t0))
    mA[(r1)*(n1)+t,(1+t*n):((t+1)*nr+1)] < 
        with(df[t+1,], c(24*t1, 6*t0))
}
    
# bond price fitting constraints
= 5for(t in 1:n) {
    mA[(r1)*(n1)+t,(1+(t1)*n):(t*n)] < 
        with(df[t+1,],c(dt5/5,dt4/4,dt3/3,dt2/2,dt1))
}
    
# additional four constraints
= (n1)*4 + n+1; c = n  ; mA[r,c] < 1
= (n1)*4 + n+2; c = n1; mA[r,c] < 1
 
= (n1)*4 + n+3; c = (n*4+1):(n*4+4)
mA[r,c] < with(df[n+1,], c(4*t3, 3*t2, 2*t1, t0))
 
= (n1)*4 + n+4; c = (n*4+1):(n*4+3)
mA[r,c] < with(df[n+1,], c(12*t2, 6*t1, 2*t0))
    
# RHS vector
vC < rep(0,n*n)
vB < c(rep(0,(n1)*(n1)), df$mln[2:(n+1)], 
        df$zrc[1],0,0,0)
    
# concatenation of matrix and vector
AA = rbind(cbind(mQ, t(mA)),
           cbind(mA, matrix(0,n*n,n*n)))
BB = c(vC, vB)
    
# solve linear system by using inverse matrix
XX = solve(AA)%*%BB
 
XX
 
cs


Running the above R code, we can obtain the following results which are the same as that of the previous post.

1
2
3
4
5
6
7
8
> df.mkt[,c(“a”,”b”,”c”,”d”,”e”)]
              a            b          c             d           e
1  3.7020077927 -3.343981336  0.8481712  1.235996e-15  0.04000000
2 -0.1354629834  0.493489440 -0.5908803  2.398419e-01  0.02500988
3  0.0080049888 -0.080382448  0.2699275 -3.340300e-01  0.16847785
4 -0.0024497295  0.045074170 -0.2946273  7.950796e-01 -0.67835432
5  0.0002985362 -0.009891144  0.1176126 -5.790532e-01  1.03931174
 
cs



Monthly Spot and Forward Curves


As noted in the previous post, in a piece-wise linear form, the spot and forward rate are as follows.

f(ti)=ait4+bit3+cit2+dit+eiy(ti)=1tiij=1(15ajΔt5j+14bjΔt4j+13cjΔt3j+12djΔt2j+ejΔtj)

Using these equations, monthly forward rate and zero curve can be obtained by using the following R code.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
# save calibrated parameters
for(i in 1:n) {
    df.mkt[i,c(“a”,“b”,“c”,“d”,“e”)] < 
        XX[((i1)*n+1):(i*n)]
}
    
# monthly forward rate and spot rate
df.mm < data.frame(
    mat = seq(0,10,1/12),y = NA, fwd = NA)
 
# which segment
df.mm$seg_no < 
    apply(df.mm, 1
          function(x) min(which(x[1]<=df.mkt$mat)) )
    
# ti^n
df.mm$t5 < df.mm$mat^5
df.mm$t4 < df.mm$mat^4
df.mm$t3 < df.mm$mat^3
df.mm$t2 < df.mm$mat^2
df.mm$t1 < df.mm$mat^1
df.mm$t0 < 1
    
nr < nrow(df.mm) # number of rows
    
# dti = ti^n-(ti-1)^n
df.mm$dt5 < c(NA,df.mm$t5[2:nr]  df.mm$t5[1:(nr1)])
df.mm$dt4 < c(NA,df.mm$t4[2:nr]  df.mm$t4[1:(nr1)])
df.mm$dt3 < c(NA,df.mm$t3[2:nr]  df.mm$t3[1:(nr1)])
df.mm$dt2 < c(NA,df.mm$t2[2:nr]  df.mm$t2[1:(nr1)])
df.mm$dt1 < c(NA,df.mm$t1[2:nr]  df.mm$t1[1:(nr1)])
    
# monthly maximum smoothness forward curve
df.mm$fwd[1< df.mkt$e[1# time 0 forward rate
df.mm$y[1]   < df.mkt$e[1# time 0 yield 
    
temp_y_sum < 0
for(i in 2:nr) {
    mat     < df.mm$mat[i]
    seg_no  < df.mm$seg_no[i] # which segment
    v_tn    < df.mm[i,c(“t4”,“t3”,“t2”,“t1”,“t0”)] 
    v_dtn   < df.mm[i,c(“dt5”,“dt4”,“dt3”,“dt2”,“dt1”)] 
    v_abcde < df.mkt[seg_no, c(“a”,“b”,“c”,“d”,“e”)]
        
    # monthly maximum smoothness forward curve
    df.mm$fwd[i] < sum(v_abcde*v_tn)
    
    # monthly yield curve
    temp_y_sum < temp_y_sum +
        sum(c(1/5,1/4,1/3,1/2,1)*v_abcde*v_dtn)
    df.mm$y[i] < (1/mat)*temp_y_sum
}
     
# Draw Graph
x11(width=8, height = 6);
plot(df.mkt$mat, df.mkt$zrc, col = “red”, cex = 1.5,
     ylim=c(0,0.12), xlab = “Maturity”
     ylab = “Interest Rate”, lwd = 10,
     main = “Monthly Maximum Smoothness Forward Rates and Spot Rates”)
text(df.mkt$mat, 
     df.mkt$zrcc(0.015,0.01,rep(0.01,3)),
     labels=c(“3M”,“1Y”,“3Y”,“5Y”,“10Y”), cex= 1.5)
            
lines(df.mm$mat, df.mm$y  , col = “blue” , lwd = 5)
lines(df.mm$mat, df.mm$fwd, col = “green”, lwd = 10)
 
legend(“bottomright”,legend=c(“Spot Curve”,
       “Forward Curve”“Input Spot Rates”),
       col=c(“blue”,“green”,“red”), pch = c(15,15,16),
       border=“white”, box.lty=0, cex=1.5)
            
cs


Drawing spot and forward curve from the maximum smoothness principle results in the following figure. These spot curve is market consistent as it fits market spot rates exactly. The forward curve shows very smooth behavior.

Maximum Smoothness Forward Rates by R code



Shortcut Method


In this problem, I guess that we can use x=A1b as shortcut method. It is because of the fact that as x is determined by Axbx=A1b, λ vector is determined by some arbitrary Q matrix except for the zero or null matrix.
Qx+ATλ=0λ=(AT)1Qx

In the above equation, if A and x are known, even though we use non-zero some arbitrary Q matrix, we can get the same x with different λ. I think this is due to the fact that our problem is exactly identified with only linear equality constraints and there is no perturbation process. The following R code consists of the linear system solution by using inverse of sub matrix mA not A.


1
2
3
4
XX2 < solve(mA)%*%vB
cbind(XX[1:(n*n)],XX2,XX[1:(n*n)]XX2)
sum(abs(XX[1:(n*n)]XX2))
 
cs


Lagrange multiplier is important concept as it is interpreted as shadow price or marginal utility in financial economics. But in the yield curve fitting context, it is of less significance. If there is no additional constraints and our purpose is to obtain only x, I think it suffice to solve Ax=b. Of course if we impose additional constraints, this shortcut will be impossible. This argument is my guess so that it is subject to correction if any inconsistency is found.

Finally, we compare forward curves from 1) R optimization, 2) Excel optimization, 3) R matrix inversion, 4) Excel matrix inversion, 5) R shortcut method. As we use Excel to illustrate the structure of sub and whole matrix or vector, coefficients can be found by using Excel solver.

The following table shows that there are negligible differences among estimated coefficients from 5 methods in the numerical perspective.




Concluding Remarks


From this post, we have implemented Adams and Deventer’s maximum smoothness forward curve. This powerful method will be used for various area, for example, interest rate trading desk since this department have a need to use exact and smooth forward curve.


To leave a comment for the author, please follow the link and comment on their blog: K & L Fintech Modeling.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)