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.
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(θP–Xt)dt+ΣdWPt
B(τ)=[11−e−τ1λτ1λ1−e−τ1λτ1λ−e−τ1λ11−e−τ2λτ2λ1−e−τ2λτ2λ−e−τ2λ⋮⋮⋮11−e−τNλτNλ1−e−τ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λ2–1λ31−e−λττ+14λ31−e−2λττ]+C[12λ2+1λ2e−λτ–14λτe−2λτ–34λ2e−2λτ–2λ31−e−λττ+58λ31−e−2λττ]+D[12λτ+1λ2e−λτ–1λ31−e−λττ]+E[3λ2e−λτ+12λτ+1λτe−λτ–3λ31−e−λττ]+F[1λ2+1λ2e−λτ–12λ2e−2λτ–3λ31−e−λττ+34λ31−e−2λττ]
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=(I–e−KPΔt)θP+e−KPΔtXt−1+η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=(I–e−KPΔt)θP , ψ1=e−KPΔt for initialization, Kalman filtering is represented as the following recursive calculations. Xt|t−1=ψ0+ψ1Xt−1|t−1Vt|t−1=ψ1Vt−1|t−1ψ‘1+Σηet|t−1=yt–B(τ)Xt–A(τ)/τevt|t−1=B(τ)Vt|t−1B(τ)‘+ΣϵXt|t=Xt|t−1+Ktet|t−1Vt|t=Vt|t−1–KtB(τ)Vt|t−1
where Kt=Vt|t−1B(τ)‘ev−1t|t−1 is the Kalman gain which reflects the uncertainty of prediction and is used as a weight for updating the time t−1 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(1−e−(λ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ΛV−1. λi is the (i,i)-th element of Λ and ωij is the (i,j)-th element of Ω=V−1ΣΣ‘(V−1)‘.
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|t−1) and its uncertainties (evt|t−1). lnLt(θ)=−NT2ln(2π)–12T∑t=1ln|evt|t−1|–12T∑t=1e‘t|t−1(evt|t−1)−1et|t−1
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 <– (1–exp(–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)–(1–exp(–la*t))/(la^3*t)+ (1–exp(–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*(1–exp(–la*t))/(la^3*t)+ 5*(1–exp(–2*la*t))/(8*la^3*t)) r4<––D*(t/(2*la)+exp(–la*t)/(la^2)–(1–exp(–la*t))/(la^3*t)) r5<––E*(3*exp(–la*t)/(la^2)+t/(2*la)+t*exp(–la*t)/(la)– 3*(1–exp(–la*t))/(la^3*t)) r6<––F*(1/(la^2)+exp(–la*t)/(la^2)–exp(–2*la*t)/(2*la^2)– 3*(1–exp(–la*t))/(la^3*t)+3*(1–exp(–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]*(1–exp(–(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)+Q # predicted cov y <– m.spot[i,] # the observed yield yimplied <– B%*%Xhat+C # the model-implied yields er <– y–yimplied # 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 <– Vhat–KG%*%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.000707, –0.003399, 0.010891, # sigma 0.032577, –0.012536, –0.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=nobs–npara) 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.
We can also plot the time-varying filtered estimates of latent factors (level, slope, curvature factors).
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.