| 37 | using namespace af; |
| 38 | |
| 39 | void simulateHestonModel(af::array &xres, af::array &vres, float T, |
| 40 | unsigned int N, unsigned int R, float mu, float kappa, |
| 41 | float vBar, float sigmaV, float rho, float x0, |
| 42 | float v0) { |
| 43 | float deltaT = T / (float)(N - 1); |
| 44 | |
| 45 | af::array x[] = {af::constant(x0, R), af::constant(0, R)}; |
| 46 | af::array v[] = {af::constant(v0, R), af::constant(0, R)}; |
| 47 | |
| 48 | float sqrtDeltaT = sqrt(deltaT); |
| 49 | |
| 50 | float sqrtOneMinusRhoSquare = sqrt(1 - rho * rho); |
| 51 | |
| 52 | float mArray[] = {rho, sqrtOneMinusRhoSquare}; |
| 53 | af::array m(2, 1, mArray); |
| 54 | |
| 55 | unsigned int tPrevious = 0, tCurrent = 0; |
| 56 | af::array zeroConstant = constant(0, R); |
| 57 | |
| 58 | for (unsigned int t = 1; t < N; t++) { |
| 59 | tPrevious = (t + 1) % 2; |
| 60 | tCurrent = t % 2; |
| 61 | |
| 62 | af::array dBt = randn(R, 2) * sqrtDeltaT; |
| 63 | af::array sqrtVLag = af::sqrt(v[tPrevious]); |
| 64 | |
| 65 | x[tCurrent] = x[tPrevious] + (mu - 0.5 * v[tPrevious]) * deltaT + |
| 66 | (sqrtVLag * dBt(span, 0)); |
| 67 | af::array vTmp = v[tPrevious] + kappa * (vBar - v[tPrevious]) * deltaT + |
| 68 | sigmaV * (sqrtVLag * matmul(dBt, m)); |
| 69 | v[tCurrent] = max(vTmp, zeroConstant); |
| 70 | } |
| 71 | |
| 72 | xres = x[tCurrent]; |
| 73 | vres = v[tCurrent]; |
| 74 | } |
| 75 | |
| 76 | int main() { |
| 77 | float T = 1; |