Robust term structure estimation in developed and emerging markets

Abstract

Despite powerful advances in interest rate curve modeling for data-rich countries in the last 30 years, comparatively little attention has been paid to the key practical problem of estimation of the term structure of interest rates for emerging markets. This may be partly due to limited data availability. However, emerging bond markets are becoming increasingly important and liquid. It is, therefore, important to be understand whether conclusions drawn from developed countries carry over to emerging markets. We estimate model parameters of fully flexible Nelson–Siegel–Svensson term structures model which has become one of the most popular term structure model among academics, practitioners, and central bankers. We investigate four sets of bond data: U.S. Treasuries, and three major emerging market government bond data-sets (Brazil, Mexico and Turkey). By including both the very dense U.S. data and the comparatively sparse emerging market data, we ensure that are results are not specific to a particular data-set. We find that gradient and direct search methods perform poorly in estimating term structures of interest rates, while global optimization methods, particularly the hybrid particle swarm optimization introduced in this paper, do well. Our results are consistent across four countries, both in- and out-of-sample, and for perturbations in prices and starting values. For academics and practitioners interested in optimization methods, this study provides clear evidence of the practical importance of choice of optimization method and validates a method that works well for the NSS model.

Appendix

1.1 Global optimization algorithms

The objective of global optimization is to find the globally best solution of (possibly non-linear) models, in the presence of multiple local optima.

1.1.1 Particle swarm optimization

PSO is a population-based metaheuristic technique. It is a fast converging algorithm, is easy to implement, and has been successfully applied to optimizing various continuous nonlinear functions in other applications.²⁴ PSO has not previously been applied to term structure estimation. PSO is essentially inspired from the social behaviour of the individuals. In a simple social setting, decision process of each individual is affected by his own experiences and other individual’s experiences. Every individual keeps a memory of their best choice, as well as the overall best choice of the population. In our PSO setting, for each observation day, dataset a set of particles search for good solutions to the NSS term structure fitting optimization problem as described in Sects. 2.2 and 2.3. Each particle is a solution of the NSS optimization problem and uses its own experience and also the experience of neighbour particles to choose how to move in the search space to find a better fitting term structure.

The PSO algorithm is initialized with a population of random candidate solutions, called particles. Each particle is assigned a random location and a random velocity, and is iteratively moved through the problem space. Every particle is attracted towards the location of the best solution achieved by the particle itself and towards the location of the best solution achieved across the whole population. At each iteration, velocity vector ($v_i$) and position (parameter) vector ($\varTheta _i$) of a particle is updated as follows:

$$\begin{aligned} v_k\leftarrow & {} v_{k-1}+U(0,\phi _1) \times (pb_{k-1}-\varTheta _{k-1})+U(0,\phi _2) \times (gb-\varTheta _{k-1}) \end{aligned}$$

(15)

$$\begin{aligned} \varTheta _k\leftarrow & {} \varTheta _{k-1}+v_{k-1} \end{aligned}$$

(16)

where the $U(0,\phi _i)$ are scalars uniformly distributed in the interval $[0,\phi _i]$, $pb_i$ is the best known position of particle i, and gb is the best known position across the entire population. The parameters $\phi _1$ and $\phi _2$ denote the magnitude of the random forces in the direction of personal best $pb_i$ and swarm best gb. The components $U(0,\phi _1) \otimes (pb_i-\varTheta _i)$ and $U(0,\phi _2) \otimes (gb-\varTheta _i)$ can be interpreted as attractive forces produced by springs of random stiffness.

PSO is an unconstrained optimization algorithm, which required us to handle the constraints in the objective function. To this end we added a penalty function that is a scalar value times the square of the violation of constraint $c_{i}$ to the objective i.e

$$\begin{aligned} \begin{aligned} D_{i}(\varTheta )&= {\left\{ \begin{array}{ll} 0 &{} \hbox {if } c_{i}(\varTheta ) \le 0 \\ C(c_{i}(\varTheta ))^2 &{} \hbox {if } c_{i}(\varTheta )>0 \end{array}\right. } \end{aligned} \end{aligned}$$

(17)

where C is large scalar value. Thus, the optimization problem becomes:

$$\begin{aligned} \min _{\varTheta } f(\varTheta ) = J(\varTheta ) + \sum \limits _{k=1}^{4} D_{k}(\varTheta ) \end{aligned}$$

(18)

There are several variants of the original PSO algorithm. Shi and Eberhart (1998) developed a variant which we call PSO-W. They modify the updating of the velocity vector to place decreasing weight on the last value, and thereby increasing weight on the individual and global bests as the iteration count increases:

$$\begin{aligned} v_k\leftarrow & {} \omega v_{k-1}+U(0,\phi _1) \times (pb_{k-1}-\varTheta _{k-1})+U(0,\phi _2) \times (gb-\varTheta _{k-1}) \end{aligned}$$

(19)

$$\begin{aligned} \varTheta _k\leftarrow & {} \varTheta _{k-1}+v_{k} \end{aligned}$$

(20)

$$\begin{aligned} \omega\leftarrow & {} \omega _p \times \omega \end{aligned}$$

(21)

where $\omega$ is termed the inertia weight, and $0<\omega _p<1$ is the decay factor. Effectively, the inertia weight preserves the initial random search direction, the magnitude of which decreases over iterations.

Another variant of PSO was proposed by Pedersen and Chipperfield (2010). This variant, named as PSO-G, disregards the personal best values and focuses on the neighborhood of the population best.

Our computational experiments show that the best results were obtained when PSO-W and PSO-G were applied sequentially. This method will be referred to as Hybrid PSO for the rest of the study.

1.1.2 Simulated annealing

Simulated annealing (SA) algorithm, proposed by Kirkpatrick et al. (1983), is a probabilistic local search algorithm that looks for the minimum of an objective function using the neighborhood information of a point in the search space. The name and inspiration of the algorithm come from annealing in metallurgy, a technique involving heating and controlled cooling of a material to increase the size of its crystals and reduce their defects. The heat causes the atoms to become unstuck from their initial positions (a local minimum of the internal energy) and wander randomly through states of higher energy; the slow cooling gives them more chances of finding configurations with lower internal energy than the initial one. By analogy with this physical process, each step of the SA algorithm replaces the current solution by a random “nearby” solution, chosen with a probability that depends both on the difference between the corresponding function values and also on a global parameter T (called the temperature), that is gradually decreased during the process.

Local search algorithms usually start with a random initial solution. A neighbour of this solution is then generated by some suitable mechanism and the change in cost is calculated. Neighbors are chosen as a network topology where the local optima are shallow or a similar topology where there are many deep local minima. If a reduction in cost is found, the current solution is replaced by the generated neighbour, otherwise the current solution is retained. The process is then repeated until no further improvement can be found in the neighbourhood of the current solution and the algorithm terminates at a local minimum. In SA, sometimes a neighbour that increases the cost is accepted, to avoid becoming trapped in a local optimum. A move that does not improve the solution is accepted with a probability that decreases with iterations. Usually, the probability is selected as $e^{-\delta /T}$ where $\delta$ is the increase in the objective function value at each iteration and T is a control parameter.

1.2 Direct search algorithms

Direct search is a class of methods for solving optimization problems that does not require any information about the gradient of the objective function. Direct search algorithms search a set of points around the current point, looking for one where the value of the objective function is lower than the value at the current point.

1.2.1 Nelder–Mead method

Nelder–Mead method is based on the idea of a “simplex”, which refers to the generalized form of a triangle in $n \ge 2$ dimensions. The algorithm is initialized with $n+1$ solutions, which form a simplex in $R^n$. The midpoint of all initial solutions excluding the worst one, is computed first. The reflection of the worst solution with respect to this midpoint is then determined as a candidate solution. Based on the parameters of the algorithm, expanded and contracted points on the line connecting the worst point and the midpoint are also determined as candidate solutions. Finally, n more candidates are generated by shrinking the simplex towards the best solution, and determining the new corners of the simplex. The best solution among all candidates then replaces the worst solution, and the algorithm is iterated until the worst solution cannot be improved.

The Nelder–Mead method (Nelder and Mead 1965; Lagarias et al. 1965) uses four scalar parameters: $\rho$ (reflection), $\chi$ (expansion), $\gamma$ (contraction), $\sigma$ (shrinkage). Let the vertices are denoted as $\varTheta ^{(1)},\varTheta ^{(2)},\ldots \varTheta ^{(n+1)}$ where n is the number of parameters to be estimated.

1.2.2 Powell’s method

Powell’s algorithm (Powell (1964)) is a generalization of a straightforward and intuitive optimization algorithm, called the “taxi-cab” algorithm. The taxi-cab algorithm has n search directions, usually the unit vectors in every dimension. The algorithm starts with an initial solution and search directions are evaluated sequentially. The best possible improvement is found in a search direction (through line minimization methods such as Golden section search), and is applied before moving on to the next search direction. The innovation of the Powell’s algorithm over the taxi-cab algorithm is the idea of updating the search directions based on the history of the search. Every time an improving solution is found, one of the search directions is replaced by the vector connecting the current solution and the improving solution, as it proved to be a good search direction.

1.3 Gradient based algorithms

Gradient based optimization algorithms iteratively search for a minimum by computing (or approximating) the gradient of the objective function. We use two popular and efficient gradient based optimization algorithms, namely Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS) and Gauss–Newton method. One of the key advantage of these algorithms compared to global and direct search algorithms is their theoretical capability of using the geometry of NSS parameters space via its function gradient information.

1.3.1 Broyden–Fletcher–Goldfarb–Shanno algorithm

The Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm (see Judd (1998)) is a quasi-Newton algorithm which is an iterative procedure with a local quadratic approximation of the function. An approximation of the Hessian of the objective function is used instead of the function itself, which decreases the complexity of the algorithm. As all variants of Newton’s methods, the idea behind them is to start at a point $\varTheta _0$ and find the quadratic polynomial $J^0(\varTheta )$ that approximates $J(\varTheta )$ to its second degree Taylor expansion evaluated at $\varTheta ^0$:

$$\begin{aligned} J^0(\varTheta )\equiv J(\varTheta _0)+\nabla J(\varTheta _0)^\prime \, (\varTheta -\varTheta _0)+\frac{1}{2}(\varTheta -\varTheta _0)^\prime \, \nabla ^2 J(\varTheta _0) (\varTheta -\varTheta _0), \end{aligned}$$

(22)

where $\nabla$ is the gradient operator, such that:

$$\begin{aligned} \nabla J(\varTheta ^0)=\left( \frac{\partial J}{\partial \beta _0}(\varTheta _0),\frac{\partial J}{\partial \beta _1}(\varTheta _0),\frac{\partial J}{\partial \beta _2}(\varTheta _0), \frac{\partial J}{\partial \beta _3}(\varTheta _0), \frac{\partial J}{\partial \tau _1}(\varTheta _0),\frac{\partial J}{\partial \tau _2}(\varTheta _0)\right) \end{aligned}$$

(23)

and $H = \nabla ^2 J(\varTheta _0)$ is the Hessian matrix of $J(\varTheta _0)$, that is the symmetric matrix containing the second-order derivatives of $J$:

$$\begin{aligned} H=\left( \frac{\partial ^2 J}{\partial \varTheta \partial \varTheta '}(\varTheta _0)\right) . \end{aligned}$$

(24)

Quasi-Newton methods specify $d_k$ (the direction vector) as:

$$\begin{aligned} d_k=-H_k^{-1}\nabla J(\varTheta _{k})^\mathsf {T}, \end{aligned}$$

(25)

The step size $\alpha _k$ is obtained by line minimization, where the direction $H_k$ is a positive definite matrix, $H_k$ may be adjusted from one iteration to the next so that $d_k$ tends to approximate the Newton direction.

1.3.2 Gauss–Newton algorithm

The idea behind the Gauss–Newton algorithm is to build a derivative matrix that contains the derivatives of the function with respect to each variable. After the first iteration, the algorithm modifies the values of the initial guess, as the optimization process goes on the variable values are updated after each iteration, until the algorithm reaches a satisfactory value (from the point of view of the operator) or the error reaches its pre-defined limit. One of the necessity of the algorithm requires the function derivative and its calculation. In the NSS estimation process, the functional form is suitable for calculation of the derivative.

1.4 Implementation details

The algorithms were coded in MATLAB. For Nelder–Mead, Powell, SA and PSO algorithms the constraints have been embedded in objective function with a penalty constant $C = 1000$, whereas they have been explicitly stated for BFGS. Stopping criterion for the algorithms has been chosen as $\epsilon =10^{-8}$. Parameter set for the Nelder–Mead algorithm is chosen as $\rho =1, \chi =2, \gamma =0.5, \sigma =0.5$. For the SA algorithm, initial value of T is chosen as 1.