# Here comes a program to compute the Nelson-Siegel approach for the

# Yield Curve. It can be modified for forward rates (straightforward

# but a little more lengthy), and to other approaches including the Svenson

# Approach and/or several Spline Interpolations.

# The data set used here are an extension of the McCulloch U.S. Treasury

# term structure data (appearing in the Handbook of Monetary Economics)

# to February 1991 as used by S-Plus.

# The results of my program are compared with those obtained from

# SPlus/Finmetrics. They agree to more than 6 digets!

NelsonSiegel =

function(Yield, Maturity)

{ # A function written by Diethelm Wuertz

# Description:

# Fit the Yield Curve by the Nelson-Siegel Method

#

# Details:

# This function finds a global solution. The start values for the

# betas are solved exactly as a function of tau using OLS.

#

# Copyright:

# Diethelm Wuertz, (c) 2004 fBonds

#

# Source:

# Partial copy from 'fBonds' from 'Rmetrics' (unpublished).

# FUNCTION:

# Find Optimal Start Solution by OLS of beta's vs. Yields:

n = length(Maturity)

gmin = 1.0e99

for (i in 1:n) {

tau = Maturity[i]

x = Maturity/tau

a = matrix(rep(NA, times = 9), nrow = 3)

a[1,1] = 1

a[1,2] = a[2,1] = mean((1-exp(-x))/x)

a[1,3] = a[3,1] = mean((1-exp(-x))/x - exp(-x))

a[2,2] = mean( ((1-exp(-x))/x)^2 )

a[2,3] = a[3,2] = mean(((1-exp(-x))/x)*((1-exp(-x))/x-exp(-x)))

a[3,3] = mean(((1-exp(-x))/x - exp(-x))^2)

b = c(

mean ( Yield ),

mean ( Yield * ((1-exp(-x))/x)),

mean ( Yield * (((1-exp(-x))/x - exp(-x)))))

beta = solve(a, b)

yfit = beta[1] + beta[2]*exp(-x) + beta[3]*x*exp(-x)

fmin = sum( (Yield-yfit)^2 )

if (fmin < gmin) {

gmin = fmin

gvec = c(beta, tau)

}

}

# Function to be optimized:

fx <- function(Maturity, x) {

x[1] + x[2] * (1-exp(-Maturity/x[4]))/(Maturity/x[4]) +

x[3] *

((1-exp(-Maturity/x[4]))/(Maturity/x[4])-exp(-Maturity/x[4]))

}

func <- function(x) { sum( (Yield - fx(Maturity, x))^2 ) }

# Optimize:

fit = nlminb(objective = func, start = gvec)

fit$start = gvec

names(fit$par) = c("beta1", "beta2", "beta3", "tau")

# Plot Curve:

yfit = fx(Maturity, gvec)

plot(Maturity, Yield, ylim = c(min(c(Yield, yfit)), max(c(Yield,

yfit))),

pch = 19, cex = 0.5, main = "Nelson-Siegel" )

lines(Maturity, yfit, col = "steelblue")

# Return Value:

fit

}

Yield = c(

0.04984, 0.05283, 0.05549, 0.05777, 0.05961, 0.06102, 0.06216, 0.06314,

0.06403,

0.06488, 0.06568, 0.06644, 0.06717, 0.06786, 0.06852, 0.06913, 0.06969,

0.07020,

0.07134, 0.07205, 0.07339, 0.07500, 0.07710, 0.07860, 0.08011, 0.08114,

0.08194,

0.08274, 0.08355, 0.08434, 0.08512, 0.08588, 0.08662, 0.08731, 0.08794,

0.08851,

0.08900, 0.08939, 0.08967, 0.08980, 0.08976, 0.08954, 0.08910, 0.08843,

0.08748,

0.08626, 0.08474, 0.08291)

Maturity = c(

0.083, 0.167, 0.250, 0.333, 0.417, 0.500, 0.583, 0.667,

0.750, 0.833,

0.917, 1.000, 1.083, 1.167, 1.250, 1.333, 1.417, 1.500,

1.750, 2.000,

2.500, 3.000, 4.000, 5.000, 6.000, 7.000, 8.000, 9.000, 10.000,

11.000,

12.000, 13.000, 14.000, 15.000, 16.000, 17.000, 18.000, 19.000, 20.000,

21.000,

22.000, 23.000, 24.000, 25.000, 26.000, 27.000, 28.000, 29.000)

NelsonSiegel(Yield, Maturity)

# $par

# beta1 beta2 beta3 tau

# 8.885514e-02 -3.671098e-02 1.398381e-11 1.025116e+00

# Compare with S-Plus - (only runs under SPlus):

# ans = term.struct(Yield, Maturity, na.rm = T, method = "ns", input =

"spot")

# c(ans$coefficients, tau = ans$tau)

# b0 b1 b2 tau

# 0.08885513 -0.03671097 5.254415e-008 1.02511

RMSE =

function(Yield, Maturity, coef)

{

sum(Yield - coef[1] +

coef[2]*(1-exp(-Maturity/coef[4]))/(Maturity/coef[4]) +

coef[3] * ( (1-exp(-Maturity/coef[4]))/(Maturity/coef[4]) -

exp(-Maturity/coef[4])) )^2

}

coef = c(8.885514e-02, -3.671098e-02, 1.398381e-11, 1.025116e+00)

RMSE(Yield, Maturity, coef)

# 1.492437

coef = c(0.08885513, -0.03671097, 5.254415e-008, 1.02511)

RMSE(Yield, Maturity, coef)

# 1.492431

DW

Thomas Steiner wrote:

>Dear Kris,

>excellent, thank you.

>

>

>

>>r = nlm(f=nelson,p=c(1,10,10,10),y=fwdcrv,steptol=1e-10,iterlim=500)

>>

>>

>

>Where is the difference to nls? I tried something like

>

>nls(formula=yc$yield~nelsonsiegel(yc$time,init), data = yc, start=init)

>

>where yc is my dataframe with time and forward rates.

>Anyway, even good starting estimates seem always to produce a singular

>gradient matrix.

>

>

>

>>Note svensson is double humped so you have an additional parameter (

>>x[5] )

>>

>>

>

>I think Svensson is 6-dimensional.

>Some interpretations of the parameters of eg Nelson-Siegel can be used

>to get better initial estimates. But this is as well mentioned in the

>paper (pdf) you linked.

>

>

>

>>Hope this helps,

>>

>>

>

>Yes, a lot. BTW: How do you calculate forward rates out of SWAP (Libor) rates?

>

>Thomas

>

>_______________________________________________

>

[hidden email] mailing list

>

https://stat.ethz.ch/mailman/listinfo/r-sig-finance>

>

>

_______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-sig-finance