[R-sig-finance] fBonds svensson fitting

Thu Feb 2 16:23:18 CET 2006

# 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
>
>_______________________________________________
>R-sig-finance at stat.math.ethz.ch mailing list
>https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>
>  
>