Читать книгу Numerical Methods in Computational Finance - Daniel J. Duffy - Страница 63

We have written some simple code to implement the explicit Euler and Richardson extrapolation methods. We present it mainly for pedagogical reasons. First, we define global arrays that hold the discrete values of the three solutions (two Euler methods on two meshes and the second-order extrapolated scheme):

// Global arrays to hold (t, y) data for meshes of size dt and dt/2 // t and y arrays on mesh dt std::vector<double> tValues; std::vector<double> yValues; // t and y arrays on mesh dt/2 std::vector<double> t2Values; std::vector<double> y2Values; // 2nd order extrapolated values on mesh dt std::vector<double> yRichValues;

The code for the schemes is:

void CurveEuler() { // db/dt = F(b), b(0) given // b(n+1) = b(n) + dt * F(b(n)) tValues = CreateMesh(L, U, NT); yValues = std::vector<double>(NT+1); yValues[0] = A; auto riccatiEquation = CreateRiccati(); double yOld; double dt = (U - L) / static_cast<double>(NT); for (std::size_t j = 1; j < tValues.size(); ++j) { yOld = yValues[j-1]; yValues[j] = yOld + dt*(*riccatiEquation)(tValues[j-1], yOld); } } void CurveRichardson() { // O(dt) -> O(dt^2) // Exercise: remove duplicate code double dt = (U - L)/ double(NT); tValues = CreateMesh(L, U, NT); yValues = std::vector<double>(NT+1); yValues[0] = A; auto riccatiEquation = CreateRiccati(); double yOld; for (std::size_t j = 1; j < tValues.size(); ++j) { yOld = yValues[j-1]; yValues[j] = yOld + dt*(*riccatiEquation)(tValues[j - 1], yOld); } // dt/2 dt *= 0.5; t2Values = CreateMesh(L, U, 2*NT+1); y2Values = std::vector<double>(2*NT+1); y2Values[0] = A; for (std::size_t j = 1; j < t2Values.size(); ++j) { yOld = y2Values[j-1]; y2Values[j] = yOld + dt*(*riccatiEquation)(t2Values[j - 1], yOld); } yRichValues = std::vector<double>(NT+1); yRichValues[0] = A; for (std::size_t j = 1; j < yRichValues.size(); ++j) { yRichValues[j] = 2.0* y2Values[2*j] - yValues[j]; } }

Some test code is:

int main() { std::cout ≪ "How many steps (need many e.g. 300000?): "; std::cin ≫ NT; CurveEuler(); CurveRichardson(); std::cout ≪ "Value at T = " ≪ U ≪ " is " ≪ std::setprecision(7) ≪ yValues[yValues.size() - 1] ≪ ", press the ANY key " ≪ std::endl; std::cout ≪ "Richardson at T = " ≪ U ≪ " is " ≪RichValues[yRichValues.size() - 1] ≪ ", press the ANY key "; // Beautiful Excel output ExcelDriver xl; xl.MakeVisible(true); xl.CreateChart(tValues, yValues, std::string("Euler for Riccati")); xl.CreateChart (tValues, yRichValues, std::string("Richardson for Riccati")); return 0; }

What is the accuracy of these schemes, and how much computer effort (size of NT and the number of function calls) is needed in order to achieve a given accuracy? To this end, we take two cases whose exact solutions have been computed using Boost C++ odeint:

 Case I: , .

 Case II: , .

In general, the library needs approximately 100 function evaluations to reach this level of accuracy. The results for explicit Euler and the extrapolated scheme for cases I and II for NT = 300, 500, 1000 and 5000 are:

 Case I: (9.2746e-6, 1.113e-5), (1.0015e-5, 1.118e-5), (1.059e-5, 1.200e-5), (1.083e-5, 1.1206e-5).

 Case II: (0.0554, 0.0559), (0.0556, 0.0559), (0.0558, 0.0559), (0.05589, 0.0559).

For a more complete discussion of Boost odeint, see Duffy (2018).

We complete our discussion of finite difference methods for (3.10). The full problem is (3.11); (3.12) is a coupled system of equations, and it can be solved in different ways. The most effective approach at this stage is to use a production third-party ODE solver in C++.