5.3.1 The Vasicek model for a bond

[ 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

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)

Multiply by an arbitrary

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 . Indeed, after integration by parts of the second order term

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

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 into overlap integrals of the form , which can be calculated analytically for a homogeneous mesh. Assuming a constant volatility in (3.5#eq.7) and drift of the form , the second last term is conveniently re-written using the coefficient as

Multiply by and write all the unknown as a function of the known values

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 when . Boundary conditions have to be justified from no-arbitrage considerations: for simplicity, the yield in (2.2.2#eq.1) is here somewhat artificially associated with the spot rate and forced to zero with the Dirichlet condition . A similar reasoning justifies the Neuman condition and is here implemented for a homogeneous grid using the second order finite difference approximation [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

The matrix is tridiagonal of the form , except the last row, where an element created by the Neuman condition has to be eliminated by hand (row minus times row ) to preserve the structure of the matrix

After substituting the value for and replacing the last equation by (5.3.1#eq.6), the linear system is solved using standard LU factorization.

Decompose the solution into a superposition of finite element roof-top functions

choosing the same elements for the test function to create as many equations as unknowns. Only the nearest neighbors contribute to the overlap integrals; these can be evaluated with a combination of the trapezoidal and the mid-point rule

where a suitable choice of the

Despite the rather sophisticated derivation, this finally yields a very elegant scheme that has been implemented in the

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
either to the
`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

`TimeTheta=0.55`

)
and a tunable integration parameter (
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.
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