previous up next SYLLABUS  Previous: 5.3 Methods for bonds  Up: 5.3 Methods for bonds  Next: 5.3.2 Extensions for derivatives


5.3.1 The Vasicek model for a bond $ \spadesuit $

[ SLIDE variational form - implicit time - FEM - code - run || VIDEO modem ( 1/ 2 ) - LAN ( 1/ 2 ) - DSL ( 1/ 2)]


The finite element method offers considerable advantages in terms of flexibility and robustness at the expense of a slightly more complicated formulation. For example, the FEM method does not suffer from the numerical instability that limits the time stepping in finite difference schemes, it converges very quickly in comparison with the Monte-Carlo method using only few driving factors and the next chapter will show how the formulation can easily be extended to deal with the American exercise style. Here we use the Vasicek equation (3.5#eq.6) only to illustrate one implementation and justify the VMARKET solution for the price of a bond $ P(r,t)$ as a function of the spot rate and time. More details about the finite element method and its formulation can be found on-line.
Since the random increments of interest rates are normally distributed in the Vasicek model, there is no real advantage to transform the problem into normalized variables. We therefore start directly from the partial differential equation (3.5#eq.6)

$\displaystyle \frac{\partial P}{\partial t} +\frac{1}{2}\sigma_s^2 \frac{\partial^2 P}{\partial r^2} +(\mu_s -\lambda\sigma_s)\frac{\partial P}{\partial r} = rP$    


Multiply by an arbitrary test function $ Q(r)$ , integrate over the domain $ \Omega=[r_-;r_+]$ where the solution is sought and formulate a variational principle

$\displaystyle \int_{r_-}^{r_+} dr Q \left[ \frac{\partial P}{\partial t} +\frac... ...tial P}{\partial r} - rP \right] = 0 \qquad \forall Q \in \mathcal{C}^1(\Omega)$ (5.3.1#eq.1)


It turns out that this variational problem is equivalent to the original equation provided that the equation is satisfied for all the test functions that are ``sufficiently general''. Galerkin's method of choosing the same functional space as the solution is an excellent starting point: here it is sufficient to keep piecewise linear function $ \mathcal{C}^1(\Omega)$ . Indeed, after integration by parts of the second order term

$\displaystyle \int_{r_-}^{r_+} dr \left[ Q\frac{\partial P}{\partial t} -\frac{... ...rtial r}\right\vert _{r_-}^{r_+} = 0 \qquad \forall Q \in \mathcal{C}^1(\Omega)$ (5.3.1#eq.2)


only first order derivatives remain, which can be evaluated using linear functions. The surface term that has been produced is used to impose so-called natural boundary conditions: the contribution from the upper boundary vanishes if $ r_+$ is large enough, since $ Q(r)\sim P(r,t)\rightarrow 0$ for $ r\rightarrow +\infty$ when $ P,Q$ are chosen from the same functional space. Essential boundary conditions will be imposed on the lower boundary to normalized the discount function where $ r_-$ is small, so that the surface term can here simply be neglected. Now discretize the time into small steps using a backward difference for the first term and a partly implicit evaluation for the rest $ P=\theta P^{t-\Delta t} + \bar{\theta} P^t$ where $ \theta=1-\bar{\theta}\in[1/2; 1]$
$\displaystyle \int_{r_-}^{r_+} dr \left[ Q\frac{P^t-P^{t-\Delta t}}{\Delta t} +... ...\frac{\partial P^{t-\Delta t}}{\partial r} - rQP^{t-\Delta t} \right) + \right.$      
$\displaystyle \left. +\bar{\theta}\left( -\frac{\sigma_s^2}{2} \frac{\partial Q... ...mbda\sigma_s)Q\frac{\partial P^t}{\partial r}\qquad - rQP^t \right) \right] = 0$      
$\displaystyle \qquad \forall Q \in \mathcal{C}^1(\Omega)$          (5.3.1#eq.3)


Discretize the interest rates, decomposing the solution into the weighted sum of finite elements roof-top functions suggested in (5.3.1#fig.1). Reassemble dependencies on $ r$ into overlap integrals of the form $ <e_i(r)\vert r^k e_j^\prime(r)>$ , which can be calculated analytically for a homogeneous mesh. Assuming a constant volatility $ \sigma_s^2\equiv e$ in (3.5#eq.7) and drift of the form $ \mu_s-\lambda\sigma_s = a(b-r) +\lambda\sigma\cos(2n\pi t/T)$ , the second last term is conveniently re-written using the coefficient $ d(t)=ab-\lambda\sigma\cos(2n\pi t/T)$ as

$\displaystyle \int_{r_-}^{r_+} dr \bar{\theta}\left[\mu_s(r)-\lambda\sigma_s\ri... ...t_{r_-}^{r_+} dr e_i r e_j^\prime -d \int_{r_-}^{r_+} dr e_i e_j^\prime \right)$    


Multiply by $ -\Delta t$ and write all the unknown $ P_j^{t-\Delta t}$ as a function of the known values $ P_j^t$
$\displaystyle \sum_{j=0}^{n-1} \left[ <e_i \vert e_j> +\theta\Delta t \left( \f... ...+ d <e_i \vert e_j^\prime> + <e_i \vert re_j> \right) \right]P_j^{t-\Delta t} =$      
$\displaystyle \;\sum_{j=0}^{n-1} \left[ <e_i \vert e_j> -\bar{\theta}\Delta t \... ... e_j^\prime> + d <e_i \vert e_j^\prime> + <e_i \vert re_j> \right) \right]P_j^t$           
$\displaystyle \forall i=0,...,n-1$          (5.3.1#eq.4)


To complete the formulation, the problem has to be supplemented with a terminal condition: from the definition of the discount function (2.2.2#eq.1) this is $ P(T-t,T)\equiv 1$ when $ t=0$ . Boundary conditions have to be justified from no-arbitrage considerations: for simplicity, the yield in $ P(0,t)=\exp[-Y(t,T)(T-t)]$ (2.2.2#eq.1) is here somewhat artificially associated with the spot rate and forced to zero with the Dirichlet condition $ P(0,t)=1, \forall t$ . A similar reasoning justifies the Neuman condition $ \partial_r P(r_+,t) \approx\partial_Y P(r_+,t)=-(T-t)P(t,T), \forall t$ and is here implemented for a homogeneous grid $ r_j=jh$ using the second order finite difference approximation $ \partial_r P(r_n)\approx[3P(r_n) -4P(r_{n-1})+P(r_{n-2})]/2h$ [1].
Because of the finite extension of finite element roof-top functions overlapping only with the nearest neighbors, the linear system of equations (5.3.1#eq.4) can be cast into

$\displaystyle \sum_{j=1}^{n-1} a_{ij} f_j^{t-\Delta t} = \sum_{j=1}^{n-1} b_{ij} f_j^t \equiv c_i, \qquad i=0,...,n-1$ (5.3.1#eq.5)


The matrix $ a_{ij}$ is tridiagonal of the form $ (l_i; d_i; r_i)$ , except the last row, where an element created by the Neuman condition $ a_{n,n-2}=1$ has to be eliminated by hand (row $ n-1$ minus $ l_{n-1}$ times row $ n$ ) to preserve the structure of the matrix

\begin{displaymath}\begin{array}{rrrl} l_{n-1} f_{n-2} &+ d_{n-1} f_{n-1} &+ r_{... ...l_{n-1}h\partial_r P(r_n) <tex2html_comment_mark>98 \end{array}\end{displaymath} (5.3.1#eq.6)



After substituting the value for $ \partial_r P(r_n)$ and replacing the last equation by (5.3.1#eq.6), the linear system is solved using standard LU factorization.
Quants: linear Galerkin finite element (FEM) discretization.
Decompose the solution $ P$ into a superposition of finite element roof-top functions

$\displaystyle P(r)=\sum_{j=0}^{n-1} P_j e_j(r),\qquad\forall Q\in\{e_i(r)\},\quad i=0,...,n-1$ (5.3.1#eq.7)


$\displaystyle e_j(x)$ $\displaystyle =$ \begin{displaymath}\left\{ \begin{array}{cl} (x-x_{j-1})/(x_j-x_{j-1}) \hspace{5... ...m}& x\in[x_j; x_{j+1}]\\ 0 & \mathrm{else} \end{array} \right.\end{displaymath} (5.3.1#eq.8)


choosing the same elements for the test function $ Q$ to create as many equations as unknowns.
Figure 5.3.1#fig.1: Linear FEM approximation illustrated with a homogeneous mesh.
\includegraphics[height=6.5cm]{figs/LinFEM.eps}
Only the nearest neighbors contribute to the overlap integrals; these can be evaluated with a combination of the trapezoidal and the mid-point rule

$\displaystyle \int_{x_i}^{x_{i+1}} f(y) dy \approx (x_{i+1}-x_i)\left[\frac{p}{... ...eft[f(x_i)+f(x_{i+1})\right] + (1-p) f\left(\frac{x_i+x_{i+1}}{2}\right)\right]$ (5.3.1#eq.9)


where a suitable choice of the tunable integration parameter reproduces piecewise constant ($ p=0$ ) or linear FEM approximations ($ p=1/3$ ) or the Crank-Nicholson method ($ p=1$ ). For a homogeneous mesh $ x=x_{i+1}-x_i=ih$ , the overlap integrals can be calculated analytically and the finite contributions yield

$\displaystyle <e_i\vert x^k e_j> \equiv \int_{x_-}^{x_+} x^k \; e_i(x) e_j(x) dx$ (5.3.1#eq.10)

\begin{displaymath}\begin{array}{cclccl} <e_i \vert e_{i\mp 1} > &=& \frac{h}{4}... ... & <e_i \vert x e_{i}^\prime > &=& \frac{h}{2}(p-1) \end{array}\end{displaymath}    


Despite the rather sophisticated derivation, this finally yields a very elegant scheme that has been implemented in the VMARKET class FEMSolution.java as 
      double twopi    = 8.*Math.atan(1.);
      double runTime  = runData.getParamValue("RunTime");
      double timeStep = runData.getParamValue("TimeStep");
      double sigma    = runData.getParamValue("Volatility");
      double theta    = runData.getParamValue("TimeTheta");
      double tune     = runData.getParamValue("TuneQuad");
      double lambda   = runData.getParamValue("MktPriceRsk");
      double t        = 1.-time/runTime;               // normalized time
      double hump     = 1.7*(t-t*t*t*t*t*t);           // volatility shaping
      double ca       = runData.getParamValue(runData.meanRevVeloNm);
      double cb       = runData.getParamValue(runData.meanRevTargNm);
      double cycles   = runData.getParamValue(runData.userDoubleNm);
      double cd       = ca*cb-lambda*sigma*Math.cos(twopi*cycles*t);
      double ce       = sigma*hump; ce=ce*ce;
                                                       //--- CONSTRUCT MATRICES
      BandMatrix a = new BandMatrix(3, f.length);      // Linear problem
      BandMatrix b = new BandMatrix(3, f.length);      //  a*fp=b*f=c
      double[]   c = new double[f.length];
                                                       // Quadrature coeff
      double h,h0,h0o,h1,h1m,h1p,h2,h2o;               //  independent of i
      double t0,t0m,t0p,t1,t1m,t1p;                    //  depending on i
      h= dx[0];
      h0o= 0.25*h*(1-tune);               h0= 0.5*h*(1+tune); 
      h1m=-0.5; h1p=-h1m;                 h1= 0.;
      h2o=-1./h;                          h2= 2./h;
      for (int i=0; i<=n; i++) {
        t0m=h*h*0.125*(1-tune)*(2*i-1);
        t0p=h*h*0.125*(1-tune)*(2*i+1);   t0=h*h*0.5*i*(tune+1);
        t1m=h*h*0.25*(-2*i-tune+1);
        t1p=h*h*0.25*(-2*i+tune+1);       t1=h*h*0.5*(tune-1);
        a.setL(i,h0o + theta   *timeStep*(t0m +0.5*ce*h2o +ca*t1m +cd*h1m));
        a.setD(i,h0  + theta   *timeStep*(t0  +0.5*ce*h2  +ca*t1  +cd*h1 ));
        a.setR(i,h0o + theta   *timeStep*(t0p +0.5*ce*h2o +ca*t1p +cd*h1p));
        b.setL(i,h0o +(theta-1)*timeStep*(t0m +0.5*ce*h2o +ca*t1m +cd*h1m));
        b.setD(i,h0  +(theta-1)*timeStep*(t0  +0.5*ce*h2  +ca*t1  +cd*h1 ));
        b.setR(i,h0o +(theta-1)*timeStep*(t0p +0.5*ce*h2o +ca*t1p +cd*h1p));
      }
      c=b.dot(f);
      double dPdy0, dPdyn, c0, cn;                     //--- BC
      a.setL(0, 0.);a.setD(0, 1.);a.setR(0, 0.);c[0]=1.;//left: Dirichlet
      double a1n= a.getD(n-1) +4.*a.getL(n-1);         // right: Neuman
      double ann= a.getR(n-1) -3.*a.getL(n-1);         //   O(h^2)
      dPdyn=-time*f[n]; cn=c[n-1]-2*dx[0]*a.getL(n-1)*dPdyn;
      a.setL(n,a1n);a.setD(n,ann);a.setR(n, 0.);c[n]=cn;
      fp=a.solve3(c);                                  //--- SOLVE
      for (int i=0; i<=n; i++) {                       //--- PLOT
        gp[i]=-Math.log(fp[i])/time;                   // yield(r)
        if (time<=timeStep) f0[i]=0.;
      }
      int i=(int)((time/runTime*n));                   // yield(t)
      f0[i]=gp[n/4];
Two band matrices a,b and a vector c are created to first assemble the linear problem (5.3.1#eq.5) using the commands of the form a.setL(j,*): they define matrix elements in row $ i$ either to the Left, Right or on the Diagonal of the matrices a,b. The right hand side vector is calculated with a product between the matrix b and the discount function f that is known from the previous time step. The solution fp is computed using LU factorization and the yield is defined from the discount function (2.2.2#eq.1) ready for plotting.
The VMARKET applet below shows the result obtained for a weakly implicit scheme ($ \theta$ or TimeTheta=0.55) and a tunable integration parameter ($ p$ or TuneQuad=1.), which is equivalent to the popular Crank-Nicholson method used by the finite differences afficinados; to be financially meaningful, the solution has of course to be independent of the numerical method.
VMARKET applet:  press Start/Stop to study the numerical properties of the finite element implementation of the Vasicek equation. Vary the implicit time parameter from explicit to implicit with TimeTheta in [0.5; 1] and tune the integration TuneQuad in [0;1] to test a piecewise constant (0) or linear (0.333) FEM approximation or even a Crank-Nicholson scheme (1).

implicit finite differences but results in the same computational cost. The additional flexibility provided by a finite elements formulation is such that the Crank-Nicholson scheme should in fact be of little more than historical interest.

SYLLABUS  Previous: 5.3 Methods for bonds  Up: 5.3 Methods for bonds  Next: 5.3.2 Extensions for derivatives