Differential Evolution (DE)
Gilli et al. (2010) applied the Differential Evolution algorithm to the Nelson-Siegel (NS) or Nelson-Siegel-Svensson (NSS) model. The authors implemented and tested this optimization and showed that it is capable of reliably solving the numerical problem of the NS or NSS model.
This DE optimization algorithm can be used from the R package, NMOF, and the authors kindly provided the R code in the appendix of their working paper. This allowed me to apply this code to some data and verify that the estimation results match those obtained using existing algorithms such as the Nelder-Mead.
Following Gilli et al. (2010), there are two notational differences from what I have been accustomed to. First, maturity is expressed in years. Second, the time decay parameter is in reciprocal form, denoted as \(\frac{1}{\lambda}\). Therefore, the NS model is represented as follow.
\[\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}\]
R code
The following R code estimates the Nelson-Siegel parameters by using the DE.
#========================================================# # Quantitative Financial Econometrics & Derivatives # ML/DL using R, Python, Tensorflow by Sang-Heon Lee # # https://shleeai.blogspot.com #--------------------------------------------------------# # Estimating the Nelson-Siegel model # using Differential Evolution (DE) #========================================================# graphics.off(); rm(list = ls()) library(NMOF) # Differential Evolution #----------------------------------------------- # Objective function #----------------------------------------------- OF <- function(param, data) { y <- data$model(param, data$tm) maxdiff <- y - data$yM # max absolute error # maxdiff <- max(abs(maxdiff)) # Sum of Squares maxdiff <- sum(maxdiff^2) if (is.na(maxdiff)) maxdiff <- 1e10 maxdiff } #======================================================= # 1. Read data #======================================================= # b1, b2, b3, lambda, RMSE for comparisons ns_reg_para_rmse1 <- c( 4.26219396, -4.08609206, -4.90893865, 3.060792, 0.04883786) ns_reg_para_rmse2 <- c( 4.97628654, -4.75365297, -6.40263059, 1.651215, 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. Estimation #======================================================= #------------------------------------------------------- # NS estimation with 1st data #------------------------------------------------------- data <- list(yM = y1, # percent tm = m/12, # in years model = NS, ww = 0.1, min = c( 0,-15,-30, 0), max = c(15, 30, 30,10) ) f_penalty <- function(mP, data) { minV <- data$min maxV <- data$max ww <- data$ww ## if larger than maxV, element in A is positiv A <- mP - as.vector(maxV) A <- A + abs(A) ## if smaller than minV, element in B is positiv B <- as.vector(minV) - mP B <- B + abs(B) # ## beta 1 + beta2 > 0 C <- ww*((mP[1L, ] + mP[2L, ]) - abs(mP[1L, ] + mP[2L, ])) A <- ww * colSums(A + B) - C A } algo <- list(nP = 100L, ## population size nG = 500L, ## number of generations F = 0.50, ## step size CR = 0.99, ## prob. of crossover min = c( 0,-15,-30, 0), max = c(15, 30, 30,10), pen = f_penalty, repair = NULL, loopOF = TRUE, ## loop over popuation? yes loopPen = FALSE, ## loop over popuation? no loopRepair = TRUE, ## loop over popuation? yes printBar = FALSE) set.seed(90) deopt <- DEopt(OF = OF, algo = algo, data = data) nsdtim_out1 <- c(deopt$xbest, sqrt(deopt$OFvalue/length(m))) #------------------------------------------------------- # NS estimation with 2nd data #------------------------------------------------------- data$yM <- y2 f_penalty <- function(mP, data) { minV <- data$min maxV <- data$max ww <- data$ww ## if larger than maxV, element in A is positiv A <- mP - as.vector(maxV) A <- A + abs(A) ## if smaller than minV, element in B is positiv B <- as.vector(minV) - mP B <- B + abs(B) # beta 1 + beta2 > 0 is not applied A <- ww * colSums(A + B) A } # new penalty without beta 1 + beta2 > 0 algo$pen <- f_penalty set.seed(90) deopt <- DEopt(OF = OF, algo = algo, data = data) nsdtim_out2 <- c(deopt$xbest, sqrt(deopt$OFvalue/length(m))) #======================================================= # 3. Results and Comparisons #======================================================= ns_reg_para_rmse1 nsdtim_out1 ns_reg_para_rmse2 nsdtim_out2 | cs |
As can be seen in the following results, It seems that parameter estimates from the DE algorithm are the same as the benchmark results.
I believe this approach is highly effective in terms of avoiding local minima and achieving fast convergence.
Reference
Gilli, M., S. Große, and E. Schumann (2010). Calibrating the Nelson–Siegel–Svensson model. COMISEF Working Paper Series No. 31. \(\blacksquare\)
No comments:
Post a Comment