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 calculate DNS shock scenarios based on IAIS (2019). We implement DNS shock generating algorithm using R code with some modification.

# IAIS DNS Shock Generating Algorithm

Risk-based Global Insurance Capital Standard is for the regulatory capital requirements against the unexpected loss of insurance companies. Especially interest rate risk is measured through the effect of interest rate shock scenarios on insurance companies’ assets and liabilities.

For this purpose International Association of Insurance Supervisors (IAIS) DNS shock generating algorithm uses Dynamic Nelson-Seigel (DNS) model to calculate five interest rate shock scenarios (mean-reversion, level-up, level-down, twist-up, twist-down). These shock scenarios are defined in the form of a baseline scenario + DNS shocks. IAIS (2019) provides the calculation methodology.
Public 2019 Field Testing Technical Specifications (https://www.iaisweb.org/page/supervisory-material/insurance-capital-standard//file/82711/public-2019-iais-field-testing-technical-specifications)
At first, given estimated DNS parameters such as $$K$$, $$\theta$$ ,$$\Sigma$$, mean reversion shock($$\text{MR Shock}$$) is calculated using $$X_0=(L_0,S_0,C_0)$$ as follows. $\text{MR Shock} = (I-e^{-K})(\theta-X_0)$ Here $$X_0=(L_0, S_0, C_0)$$ is estimated using the Nelson-Siegel model with the year-end term structure of interest rates using nonlinear least squares. But since the DNS model is typically estimated using the Kalman filter, we can use the last latent factor estimates. $$I$$ is an identity matrix.

To define the remaining four shocks, we introduce $$M$$ as follows. \begin{align} M =\sqrt{ (\Sigma \Sigma^T) \odot \left( \frac{1-e^{-(K_i + K_j)}}{K_i + K_j} \right)_{ij} } \end{align} Next, we define $$N$$ for selecting the first and second major factors according to the magnitude of eigenvalue from principal component analysis (PCA). \begin{align} N = \begin{bmatrix} LOT&0&0\\ 0&a&0\\ 0&0&b \end{bmatrix} \end{align} Here, $$a=\sum_{\tau=1}^{\tau=LOT} \left[ \frac{1-e^{-\lambda\tau}}{\lambda\tau}\right]$$, $$b=\sum_{\tau=1}^{\tau=LOT} \left[ \frac{1-e^{-\lambda\tau}}{\lambda\tau} -e^{-\lambda\tau} \right]$$, and $$LOT$$ denotes long-term maturity (i.e. 20-year).

If we apply eigen decomposition to diagonalize $$N^{‘}N$$, we can get two eigenvectors, $$e_1$$ and $$e_2$$, which are corresponding to the largest and the second largest eigenvalues respectively. $$e_3$$ is discarded becuase its eigenvalue is the smallest. $$e_1$$ and $$e_2$$ are orthonormal, which means $$\lVert e_1 \rVert = \lVert e_2 \rVert =1$$ with the Euclidean distance of 1.

$$Me_1$$ and $$Me_2$$ represent level and twist shock respectively. To make the sum of the term structure of yields only from the twist shock to zero, we apply the affine rotation to $$Me_1$$ and $$Me_2$$ simultaneously and redefine the level and twist shock as follows.

\begin{align} Twist &= cos(\theta)Me_2 – sin(\theta)Me_1, \\ Level &= cos(\theta)Me_1 + sin(\theta)Me_2 \end{align} If $$S_1(\tau)$$ and $$S_2(\tau)$$ are the term structure of shock from $$Me_1$$ and $$Me_2$$ across maturities($$\tau$$), the following relationship is established. \begin{align} &\sum_{\tau=1}^{LOT} (cos(\theta)S_2(\tau) – sin(\theta)S_1(\tau)) = 0 \\ &\rightarrow tan(\theta) = \frac{\sum_{\tau=1}^{LOT} S_2(\tau)}{\sum_{\tau=1}^{LOT} S_1(\tau)} \end{align} Solving for $$\theta$$, the twist and level shock at the 99.5% are calculated with the following expressions. \begin{align} \text{Twist}^{99.5th} &= N^{-1}(0.995) \times Twist, \\ \text{Level}^{99.5th} &= N^{-1}(0.995) \times Level \end{align} Finally, we can obtain the following five DNS shocked curves by adding the above shocks to the initial yield curve (baseline). \begin{align} \text{Mean-Reversion} &= \text{Base} + X_0 \times \text{MR Shock}, \\ \text{Level-Up(Down)} &= \text{Base} \pm X_0 \times \text{Level}^{99.5th}, \\ \text{Twist-Up(Down)} &= \text{Base} \pm X_0 \times \text{Twist}^{99.5th} \end{align} Here $$X_0=(L_0,S_0,C_0)$$ and $$B(\tau)$$ is the factor loading matrix of DNS model and $$Base$$ denotes initial spot yield curve.

This is all about IAIS(2019) DNS shock generating algorithm but we think that the above equation for $$M$$ is needed to be more accurate.

The mathematical inconsistency of the above expression can cause confusion for practitioners in insurance companies. When we calculate an exponential of a matrix, we should use the matrix exponential function, not scalar exponential one. In R, we should expm() not exp(). Also, it is not typical to use the expression of a square root of a matrix. In this case, we’d rather use the Cholesky decomposition of a matrix than a square root of it. Therefore, it is necessary to represent M in a mathematically consistent way as follows.
\begin{align} &M = \text{lower triangular matrix of } \\ &\textbf{Chol}\left( (\Sigma \Sigma^T) \odot \left[ \frac{1-{(e^{-K})}_i \times {(e^{-K})}_j}{K_i + K_j} \right]_{ij} \right) \end{align}

The following R code implements the IAIS DNS shock generating algorithm with given DNS parameter estimates.

 #=========================================================================## Financial Econometrics & Derivatives, ML/DL using R, Python, Tensorflow # by Sang-Heon Lee## https://kiandlee.blogspot.com#————————————————————————-## IAIS DNS Shock Generating Algorithm#=========================================================================#        library(expm) # for matrix exponential        graphics.off()  # clear all graphs    rm(list = ls()) # remove all files from your workspace        # DNS model parameter estimates using Kalman filter (2011/07~2019/06)            nk <– 3  # number of factor                # mean-resversion speed parameters        Kappa = diag(c(0.09916472,0.571726962,0.637339973))                # long-term mean        Theta = c(0.036545095,–0.014733322,–0.007498777)                # lower triangular volatility matrix        Sigma = rbind(c(0.005454287 ,  0          ,           0),                      c(–0.004662791,  0.002776609,           0),                      c(–2.4875E–05 , –0.001050356, 0.009006942))                 lambda=0.365916203              # time decay parameter        X0=c(0.02024,–0.00420,–0.00791) # last factor estimates        # BB = DNS factor loading matrix        LOT    <– 20        v.1LOT <– 1:LOT # 1,2,3,…,20        BB     <– t(rbind(rep(1,LOT),                  (1–exp(–lambda*v.1LOT))/(lambda*v.1LOT),                 (1–exp(–lambda*v.1LOT))/(lambda*v.1LOT)–                  exp(–lambda*v.1LOT)))        # Definitions of M and N, calculation of  N’N                 eK <– expm(–Kappa)        m.temp <– matrix(0,nk,nk)        for (i in 1:nk) {            for (j in 1:nk) {                m.temp[i,j] <– (1–eK[i,i]*eK[j,j])/                              (Kappa[i,i]+Kappa[j,j])            }        }                    M  <– t(chol((Sigma%*%t(Sigma))*m.temp))        N  <– diag(c(LOT,sum(BB[,2]), sum(BB[,3])))%*%M        NN <– t(N)%*%N        # Eigen Decomposition        eigNN <– eigen(NN, symmetric = TRUE)        v     <– eigNN$values e <– eigNN$vectors        # Me1, Me2        Me1 <– M%*%e[,1]; Me2 <– M%*%e[,2]        # S1(1:LOT), S2(1:LOT)         S1 <– BB%*%Me1; S2 <– BB%*%Me2            # calculation of angle        angle <– atan(sum(S2)/sum(S1))            # Mean-Reversion Shock        MRS.shock   <– (diag(nk)–expm(–Kappa))%*%(Theta–X0)        # Level and Twist Shock        NI95  <– qnorm(0.995,0,1)        Level.shock <– NI95*(cos(angle)*Me1+sin(angle)*Me2)        Twist.shock <– NI95*(cos(angle)*Me2–sin(angle)*Me1)        # market spot rates and maturities        tau  <– c(0.25, 0.5, 0.75, 1, 1.5, 2,                  2.5, 3, 4, 5, 7, 10, 20)                spot <– c(0.015241, 0.016393, 0.017965,                   0.018897, 0.020274, 0.021070,                   0.021723, 0.021813, 0.023859,                   0.024832, 0.025135, 0.024984, 0.025005)            # mB = DNS factor loading matrix with market maturity        mB <– t(rbind(rep(1,length(tau)),               (1–exp(–lambda*tau))/(lambda*tau),               (1–exp(–lambda*tau))/(lambda*tau)–exp(–lambda*tau)))            # DNS shock size with market maturities        MRS.add   <– mB%*%MRS.shock        Level.add <– mB%*%Level.shock        Twist.add <– mB%*%Twist.shock            # DNS Shocked Curve = Initial Spot Curve + DNS shock        MRS.curve      <– spot + MRS.add        Level.Up <– spot + Level.add        Level.Dn <– spot – Level.add        Twist.Up <– spot + Twist.add        Twist.Dn <– spot – Twist.add            # DNS Shock graph        x11(width=6, height=5);        plot(x=tau, y=MRS.add,type=“b”,col=“darkgreen”,lty=1,             ylab=“DNS Shock”,lwd=2,xlab=“Maturity”,             main=“IAIS DNS Shocks”,             ylim=c(min(–Level.add),max(Level.add)))                lines(x=tau, y= Level.add, type=“b”,col=“red” ,lty=1,lwd=2)        lines(x=tau, y=–Level.add, type=“b”,col=“blue”,lty=1,lwd=2)        lines(x=tau, y= Twist.add, type=“b”,col=“red” ,lty=1,lwd=2,pch=16)        lines(x=tau, y=–Twist.add, type=“b”,col=“blue”,lty=1,lwd=2,pch=16)                legend(“topright”, pch=c(1,1,1,16,16), cex = 0.9,          col   =c(“darkgreen”,“red”,“blue”,“red”,“blue”),          legend=c(“MRS”,“Level Up”,“Level Down”,“Twist Down”,“Twist Up”))        # DNS Shocked Curve graph        x11(width=6, height=5);        plot(x=tau, y=spot,type=“b”,col=“black”,lty=1,             ylab=“DNS Shocked Curve”, lwd=2,xlab=“Maturity”,             main=“IAIS DNS Shocked Curves”,             ylim=c(min(Level.Dn),max(Level.Up)), pch=16)                lines(x=tau, y=MRS.curve ,type=“b”,col=“darkgreen”,lty=1,lwd=2)        lines(x=tau, y=Level.Up, type=“b”,col=“red” ,lty=1,lwd=2)        lines(x=tau, y=Level.Dn, type=“b”,col=“blue”,lty=1,lwd=2)        lines(x=tau, y=Twist.Up, type=“b”,col=“red” ,lty=1,lwd=2,pch=16)        lines(x=tau, y=Twist.Dn, type=“b”,col=“blue”,lty=1,lwd=2,pch=16)                legend(“topright”,pch=c(16,1,1,1,16,16), cex = 0.9,               col   =c(“black”,“darkgreen”,“red”,“blue”,“red”,“blue”),               legend=c(“Baseline”,“MRS”,“Level Up”,                        “Level Down”,“Twist Down”,“Twist Up”))        Colored by Color Scripter cs

Running this R code produces the two figures. the first figure draws five DNS shocks.

The second figure draws five DNS shocked curves (scenarios) and a baseline curve.

This completes the implementation of IAIS DNS shock generating algorithm using R.

The remaining job is two-fold. The one is to apply the Smith-Wilson extrapolation to these DNS shock scenarios. The other is to apply the Hull-White 1 factor simulation model with the output of the former as the initial yield curve. the output of the former refers to as the deterministic DNS shock scenario and the latter the stochastic DNS shock scenario. These two works using R are explained in our next posts.  $$\blacksquare$$