Arbitrage-free Nelson Siegel model with R code

[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.


Sang-Heon Lee

This article explains how to estimate parameters of the Arbitrage-Free dynamic Nelson-Siegel (AFNS) model (Christensen, Diebold, and Rudebusch; 2009, Christensen, Diebold, and Rudebusch; 2011) using Kalman filter. We estimate not only parameters but also filtered latent factor estimates such as level, slope, and curvature using R code.




1. AFNS model


AFNS model can be expressed as a state state model which consists of measurement equation and state equation as follows.
yt(τ)=B(τ)Xt+A(τ)τ+ϵt(τ)dXt=KP(θPXt)dt+ΣdWPt

B(τ)=[11eτ1λτ1λ1eτ1λτ1λeτ1λ11eτ2λτ2λ1eτ2λτ2λeτ2λ11eτNλτNλ1eτNλτNλeτNλ],
yt(τ)=[yt(τ1)yt(τ2)yt(τN)],ϵt(τ)=[ϵt(τ1)ϵt(τ2)ϵt(τN)],
KP=[κP11000κP22000κP33],θP=[θP1θP2θP3],Σ=[σ1100σ21σ220σ31σ32σ33],dWPt=[dWP1tdWP2tdWP3t]
Christensen et al. (2009, 2011) introduces the yield adjustment term to ensure the no-arbitrage condition for bond price. This yield adjustment term looks very complicated but can be calculated using estimated parameters such as τ, λ, Σ (see Christensen et al. (2009, 2011) for the detailed derivation).
A(τ)τ=Aτ26+B[12λ21λ31eλττ+14λ31e2λττ]+C[12λ2+1λ2eλτ14λτe2λτ34λ2e2λτ2λ31eλττ+58λ31e2λττ]+D[12λτ+1λ2eλτ1λ31eλττ]+E[3λ2eλτ+12λτ+1λτeλτ3λ31eλττ]+F[1λ2+1λ2eλτ12λ2e2λτ3λ31eλττ+34λ31e2λττ]
where A=coef(1,1),B=coef(2,2),C=coef(3,3),D=coef(1,2),E=coef(1,3),F=coef(2,3) if we denote coef(i,j)=σi1σj1+σi2σj2+σi3σj3.

Since AFNS model have a continuous-time representation, we should express AFNS model as the discretized form for estimating parameters using discrete time series data as follows. yt(τ)=B(τ)Xt+A(τ)τ+ϵt(τ)Xt=(IeKPΔt)θP+eKPΔtXt1+ηt
where Δt is 1/12 for monthly data and 1/52 for weekly data.

AFNS model as a state space model is linear with respect to factors, we can use Kalman filter to estimate parameters. If we use ψ0=(IeKPΔt)θP , ψ1=eKPΔt for initialization, Kalman filtering is represented as the following recursive calculations. Xt|t1=ψ0+ψ1Xt1|t1Vt|t1=ψ1Vt1|t1ψ1+Σηet|t1=ytB(τ)XtA(τ)/τevt|t1=B(τ)Vt|t1B(τ)+ΣϵXt|t=Xt|t1+Ktet|t1Vt|t=Vt|t1KtB(τ)Vt|t1
where Kt=Vt|t1B(τ)ev1t|t1 is the Kalman gain which reflects the uncertainty of prediction and is used as a weight for updating the time t1 prediction after observing time t data.

The conditional covariance (Ση) and unconditional covariance (Ση,0) are defined as Ση=Δt0exp(Kps)ΣΣexp(Kps)dsΣη,0=0exp(Kps)ΣΣexp(Kps)ds
Christensen and Rudebusch (2015) explain the way to calculate the analytical covariance as follows. Ση=V(ωijλi+λj(1e(λi+λj)Δt))n×nV
where V and Λ are the eigenvector matrix and the diagonal matrix with eigenvalues respectively from the diagonalization of KP=VΛV1. λi is the (i,i)-th element of Λ and ωij is the (i,j)-th element of Ω=V1ΣΣ(V1).

From the equation for the conditional covariance matrix, we can eaily derive the unconditional covariance matrix. Ση,0=V(ωijλi+λj)n×nV
Since Kalman fitler is an iterative procedure, we use X1|0=θP and V1|0=Ση,0 for initial guess. We use numerical optimization algorithm to search parameters for maximizing the log likelihood function that is constructed from conditional prediction errors (et|t1) and its uncertainties (evt|t1). lnLt(θ)=NT2ln(2π)12Tt=1ln|evt|t1|12Tt=1et|t1(evt|t1)1et|t1




2. R code for AFNS model


Now we can implement the AFNS model using R code as follows. Our estimation is based on monthly KTB (Korean Treasury Bond) from January 2011 to December 2019. The data are end-of-month, zero-coupon yields at thirteen maturities: 0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 4, 5, 7, 10, and 20 years.


#=========================================================================#
# Financial Econometrics & Engineering, ML/DL using R, Python, Tensorflow 
# by Sang-Heon Lee 
#
# https://kiandlee.blogspot.com
#————————————————————————-#
# AFNS estimation using Kalman filter
#=========================================================================#
 
library(readxl)
library(expm)
library(numDeriv)
 
graphics.off()  # clear all graphs
rm(list = ls()) # remove all files from your workspace
 
setwd(“D:/a_book_FIER_Ki_Lee/ch04_KICS/code/dns_afns_est”)
 
# DNS factor loading matrix
NS.B<function(lambda, tau)
{
    col1 < rep.int(1,length(tau))
    col2 < (1exp(lambda*tau))/(lambda*tau)
    col3 < col2  exp(lambda*tau) 
    return(cbind(col1,col2,col3))
}
 
# yield adjustment term in AFNS
AFNS.C<function(sigma, lambda, tau)
{
    s < sigma; la < lambda; t < tau
    
    s11<s[1,1]; s12<s[1,2]; s13<s[1,3]
    s21<s[2,1]; s22<s[2,2]; s23<s[2,3]
    s31<s[3,1]; s32<s[3,2]; s33<s[3,3]
  
    A<s11^2+s12^2+s13^2; D<s11*s21+s12*s22+s13*s23
    B<s21^2+s22^2+s23^2; E<s11*s31+s12*s32+s13*s33 
    C<s31^2+s32^2+s33^2; F<s21*s31+s22*s32+s23*s33
    
    r1<A*t^2/6
    r2<B*(1/(2*la^2)(1exp(la*t))/(la^3*t)+
                       (1exp(2*la*t))/(4*la^3*t))
    r3<C*(1/(2*la^2)+exp(la*t)/(la^2)t*exp(2*la*t)/(4*la)
           3*exp(2*la*t)/(4*la^2)2*(1exp(la*t))/(la^3*t)+
           5*(1exp(2*la*t))/(8*la^3*t))
    r4<D*(t/(2*la)+exp(la*t)/(la^2)(1exp(la*t))/(la^3*t))
    r5<E*(3*exp(la*t)/(la^2)+t/(2*la)+t*exp(la*t)/(la)
           3*(1exp(la*t))/(la^3*t))
    r6<F*(1/(la^2)+exp(la*t)/(la^2)exp(2*la*t)/(2*la^2)
           3*(1exp(la*t))/(la^3*t)+3*(1exp(2*la*t))/(4*la^3*t))
    return(r1+r2+r3+r4+r5+r6)
}
 
# parameter restrictions
trans<function(b)
{
    bb < b
    bb[1]  < 1/(1+exp(b[1]))  # kappa11
    bb[13< b[13]^2          # lambda
    bb[14:npara] < b[14:npara]^2          # measurement error
    return(bb)
}
 
# log likelihood function
loglike<function(para_un,m.spot)
{
    # parameter restrictions
    para < trans(para_un)
  
    # restricted parameters
    kappa  < rbind(c(para[1],0,0),
                    c(0,para[2],0),
                    c(0,0,para[3]))
    sigma  < rbind(c(para[4],0,0),
                    c(para[5],para[6],0),
                    c(para[7],para[8],para[9]))
    theta  < para[10:12]
    lambda < para[13]
    H      < diag(para[14:npara])
   
    B  < NS.B(lambda,v.mat); tB < t(B) # factor loading matrix
    C  < AFNS.C(sigma,lambda,v.mat)     # yield adjustment
  
    # Conditional and Unconditional covariance matrix : Q, Q0
    m    < eigen(kappa) 
    eval < m$values 
    evec < m$vectors; ievec<solve(evec)
    Smat < ievec%*%sigma%*%t(sigma)%*%t(ievec)
    Vdt  < Vinf < matrix(0,nf,nf)
    
    for(i in 1:nf) { for(j in 1:nf) {
        Vdt[i,j] <Smat[i,j]*(1exp((eval[i]+eval[j])*dt))/
                             (eval[i]+eval[j]) # conditional
        Vinf[i,j]<Smat[i,j]/(eval[i]+eval[j]) # unconditional
    }}
    
    # Analytical Covariance matrix
    # Q : conditional, Q0 : unconditional
    Q  < evec%*%Vdt%*%t(evec)
    Q0 < evec%*%Vinf%*%t(evec)
    
    # initialzation of vector and matrix
    prevX < theta; prevV < Q0
    Phi1  < expm(kappa*dt)
    Phi0  < (diag(nf)Phi1)%*%theta
    loglike < 0 # log likelihood function
    
    for(i in 1:nobs) {
        Xhat < Phi0+Phi1%*%prevX        # predicted state
        Vhat < Phi1%*%prevV%*%t(Phi1)+# predicted cov
        
        y        < m.spot[i,] # the observed yield
        yimplied < B%*%Xhat+# the model-implied yields
        er       < yyimplied # prediction error
 
        # updating 
        ev < B%*%Vhat%*%tB+H; iev<solve(ev)
        KG < Vhat%*%tB%*%iev # Kalman Gain
        
        prevX < Xhat+KG%*%er       # E[X|y_t]   updated state 
        prevV < VhatKG%*%B%*%Vhat # Cov[X|y_t] updated cov
        
        # log likelihood function
        loglike < loglike  0.5*(nmat)*log(2*pi)
                             0.5*log(det(ev))0.5*t(er)%*%iev%*%er
        gm.factor[i,] << prevX
    }
    
    return(loglike)
}
 
 
#=========================================================================#
#  Main : AFNS term structure model estimation
#=========================================================================#
 
    dt < 1/12 #1/52 # weekly
    nf < 3
    
    # read excel spot data
    file < “spot_2011_2019.xlsx”; sheet < “monthly”
    df.spot < read_excel(file,sheet,col_names=TRUE)
    
    # divide date and data
    v.ymd  < df.spot[,1]
    v.mat  < as.numeric(colnames(df.spot)[1])
    m.spot < as.matrix(df.spot[,1])
    
    nmat < length(v.mat) # # of maturities
    nobs < nrow(m.spot)  # # of observations
    
    # factor estimates
    gm.factor < matrix(0,nobs,nf)
    
    #—————————————————–
    # initial guess for unconstrained parameters(para_un)
    #—————————————————–
    init_para_un < c(
      1.226637,  0.840692,  0.603496# kappa
      0.006327,  
     0.005464,  0.003441,  
     0.0007070.003399,  0.010891# sigma
      0.0325770.0125360.002748# theta 
      0.5    ,                       # lambda
      rep(0.000524,nmat)              # measurement error
    )
    
    npara < length(init_para_un) # # of observations
    
    m<optim(init_para_un,loglike,
             control = list(maxit=5000, trace=2),
             method=c(“Nelder-Mead”),m.spot=m.spot)
    prev_likev < m$value
    v.likev    < m$value
    
    i < 1
    repeat {
        print(paste(i,“-th iteration”))
        m<optim(m$par,loglike,
                 control = list(maxit=5000, trace=2),
                 method=c(“Nelder-Mead”),m.spot=m.spot)
        
        v.likev < cbind(v.likev, m$value)
        print(paste(m$value,” <- likelihood value"))
        
        if (abs(m$value  prev_likev) < 0.00000000001) {
            break
        }
        prev_likev < m$value
        i < i + 1
        
        print(v.likev)
    }
    
    name_theta < 
      c(“kappa_1”   , “kappa_2”   , “kappa_3” ,
        “sigma_11”  , “sigma_21”  , “sigma_22”,
        “sigma_31”  , “sigma_32”  , “sigma_33”,
        “theta_1”   , “theta_2”   , “theta_3” ,
        “lambda”    , “epsilon_1” ,
        “epsilon_2” , “epsilon_3” , “epsilon_4”,
        “epsilon_5” , “epsilon_6” , “epsilon_7”,
        “epsilon_8” , “epsilon_9” , “epsilon_10”,
        “epsilon_11”“epsilon_12”“epsilon_13”)
 
    # draw 3 factor estimates
    x11(width=6, height=5);
    matplot(gm.factor,type=“l”, ylab=“L,S,C”, lty = 1
            main = “AFNS 3 Factor Estimates (L,S,C)”, lwd=2)
 
    # Delta method for statistical inference
    grad    < jacobian(trans, m$par)
    hess    < hessian(func=loglike, x=m$par,m.spot=m.spot)
    vcv_con < grad%*%solve(hess)%*%t(grad)
    
    # parameter | std.err | t-value | p-value
    theta   < trans(m$par)
    stderr  < sqrt(diag(vcv_con))
    tstat   < theta/stderr
    pvalue  < 2*pt(abs(tstat),df=nobsnpara)
    df.est  < cbind(theta, round(stderr,4),
                     round(tstat,4), round(pvalue,4))
 
    rownames(df.est) < name_theta # parameter name
    colnames(df.est) <c(“parameter”,“std.err”,“t-stat”,“p-value”)
    print(df.est)
cs


3. Estimated AFNS model


Running the above R code for the AFNS model, we can get the estimated parameters and the latent factor estimates(L,S,C). The following figure prints out the convergence of the log-likelihood function and estimated parameters with standard errors, t-statistics, and p-values. We use the delta method for statistical inference.

AFNS model estimation r code

We can also plot the time-varying filtered estimates of latent factors (level, slope, curvature factors).

level slope curvature r code afns

Reference

Christensen, J. H. E., F. X. Diebold, and G. D. Rudebusch, 2009, An Arbitrage-free Generalized Nelson-Siegel Term Structure Model, Econometrics Journal 12, 33~64.

Christensen, J. H. E., F. X. Diebold, and G. D. Rudebusch, 2011, The Arbitrage-free Class of Nelson-Siegel Term Structure Models. Journal of Econometrics 164, 4-20.

Christensen, J. H. E. and G. D. Rudebusch (2015), Analytical Formulas for the Second Moment in Affine Models with Stochastic Volatility, working paper


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.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)