Site icon R-bloggers

Non-linear Optimization by using constrOptim.nl R function

sang-heon lee
[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.
< !--shlee --> This post shows how to use constrOptim.nl() R function to solve non-linear optimization problem with or without equality or inequality constraints. Nelson-Siegel yield curve model is used as an target example. Its calculation time is faster than nloptr() function.
< !--
–>

Nelson-Siegel model using constrOptim.nl() R function



In this post, the non-linear least squares problem is solved by using constrOptim.nl() function from alabama R package. Like the previous posts, Nelson-Siegel model is also used for a testing purpose.

For your information, this blog has dealt with Nelson-Siegel model as follows.




Nelson-Siegel model as a non-linear optimization problem


Nelson-Siegel model is a non-linear least square problem with 4 parameters with some inequality constraints.
\[\begin{align} y(\tau) &= \beta_1 + \beta_2 \left( \frac{1-e^{- \tau \lambda }}{\tau \lambda }\right) \\ &+ \beta_3 \left(\frac{1-e^{- \tau \lambda }}{\tau \lambda }-e^{- \tau \lambda }\right) \end{align}\] \[\begin{align} \beta_1 > 0, \beta_1 + \beta_2 >0, \lambda > 0 \end{align}\]
Here, \(\tau\) is a maturity and \(y(\tau)\) is a continuously compounded spot rate with \(\tau\) maturity. \(\beta_1, \beta_2, \beta_3\) are coefficient parameters and \(\lambda\) is the time decay parameter.



constrOptim.nl() function


After installing alabama R package, constrOptim.nl() function can be used. Its arguments are so straightforward to understand that I just point out two arguments for constraints. The following two arguments denote inequality and equality constraint functions respectively. In particular, unlike nloptr() function, the inequality constraint function (hin) returns a vector of greater than zero constraints. You can easily find their functional form in the R code which will be shown below


# hin : a vector function specifying inequality constraints such that hin[j] > 0 for all j

# heq : a vector function specifying equality constraints such that heq[j] = 0 for all j



For further information about its settings, please see the document (https://www.rdocumentation.org/packages/alabama/versions/2022.4-1/topics/constrOptim.nl).

< !--콘텐츠 내 자동 삽입 광고 배치하기 (compact)-->< !--–>

R code


The following R code is similar to the previous post in that it estimates parameters of Nelson-Siegel model with yield curves at two dates.

< !--shlee -->
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
122
123
#========================================================#
# Quantitative ALM, Financial Econometrics & Derivatives 
# ML/DL using R, Python, Tensorflow by Sang-Heon Lee 
#
# https://kiandlee.blogspot.com
#——————————————————–#
# Nelson-Siegel using constrOptim.nl function
#========================================================#
 
graphics.off(); rm(list = ls())
 
library(alabama)
    
#———————————————–
# constrOptim.nl objective function
#———————————————–
    objfunc_co <– function(para, y, m) {
        beta <– para[1:3]; la <– para[4]
        C  <– cbind(
            rep(1,length(m)),    
            (1–exp(–la*m))/(la*m),             
            (1–exp(–la*m))/(la*m)–exp(–la*m))
        return(sqrt(mean((y – C%*%beta)^2)))
    }
    
#———————————————–
# constrOptim.nl constraint function ( > : larger )
#———————————————–
    # beta0 + beta1 > 0
    constfunc_co <– function(para, y, m) {
        beta  <– para[1:3]; la <– para[4]
        
        h <– rep(NA, 4)
        h[1] <– beta[1]         # b1 > 0
        h[2] <– beta[1]+beta[2] # b1+b2 > 0
        h[3] <– la–0.001        # lambda > 0.001
        h[4] <– 0.1–la          # lambda < 0.1
                   
        return(h)
    }
 
#=======================================================
# 1. Read data
#=======================================================
    
    # b1, b2, b3, lambda for comparisons
    ns_reg_para_rmse1 <– c(
        4.26219396, –4.08609206, –4.90893865,  
        0.02722607,  0.04883786)
    ns_reg_para_rmse2 <– c(
        4.97628654, –4.75365297, –6.40263059,  
        0.05046789,  0.04157326)
    
    str.zero <– “
        mat     rate1      rate2
        3       0.0781     0.0591
        6       0.1192     0.0931
        9       0.1579     0.1270
        12      0.1893     0.1654
        24      0.2669     0.3919
        36      0.3831     0.8192
        48      0.5489     1.3242
        60      0.7371     1.7623
        72      0.9523     2.1495
        84      1.1936     2.4994
        96      1.4275     2.7740
        108     1.6424     2.9798
        120     1.8326     3.1662
        144     2.1715     3.4829
        180     2.5489     3.7827
        240     2.8093     3.9696″
    
    df <– read.table(text = str.zero, header=TRUE)
    m  <– df$mat
    y1 <– df$rate1; y2 <– df$rate2
 
#==============================================================
# 2. constrOptim.nl 
#  : Nonlinear least squares with constraints
#==============================================================
    
    #———————————————
    # NS estimation with 1st data
    #———————————————
    y <– y1
    x_init <– c(y[16], y[1]–y[16], 
                2*y[6]–y[1]–y[16], 0.0609)
    
    m_coptim <– constrOptim.nl(par = x_init, y = y, m = m,
        fn = objfunc_co, hin = constfunc_co,
        control.optim = list(maxit=50000, trace=0, 
            reltol = 1e–16, abstol = 1e–16) ,
        control.outer = list(eps = 1e–16, itmax = 50000, 
            method = “Nelder-Mead”, NMinit = TRUE))
    
    ns_coptim_out1 <– c(m_coptim$par, m_coptim$value)
    
    #———————————————
    # NS estimation with 2nd data
    #———————————————
    y <– y2
    x_init <– c(y[16], y[1]–y[16], 
                2*y[6]–y[1]–y[16], 0.0609)
    
    m_coptim <– constrOptim.nl(par = x_init, y = y, m = m,
        fn = objfunc_co, hin = constfunc_co,
        control.optim = list(maxit=50000, trace=0, 
            reltol = 1e–16, abstol = 1e–16) ,
        control.outer = list(eps = 1e–16, itmax = 50000, 
            method = “Nelder-Mead”, NMinit = TRUE))
    
    ns_coptim_out2 <– c(m_coptim$par, m_coptim$value)
    
#=======================================================
# 3. Results and Comparisons
#=======================================================
    
    ns_reg_para_rmse1
    ns_coptim_out1
    
    ns_reg_para_rmse2
    ns_coptim_out2
    
cs


< !--shlee --> The following estimation results show that non-linear optimization results from constrOptim.nl() function are nearly similar to the answers.



Using these estimated parameters at two selected dates, two fitted term structure of interest rates can be generated as follows.



Concluding Remarks


This post shows how to use constrOptim.nl() R function to perform non-linear optimization problem with or without equality or inequality constraints. The reason why I introduce this optimization function is that its computation time tends to be a little faster than nloptr() function. \(\blacksquare\)


< !--shlee -->
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.