financial/heston_model.cpp
/**********************************************************************************************
* Copyright (c) 2015, Michael Nowotny
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification,
* are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation and/or other
* materials provided with the distribution.
*
* 3. Neither the name of the copyright holder nor the names of its contributors may be used
* to endorse or promote products derived from this software without specific
* prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
* A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
* OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
* TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
* PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
* NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
***********************************************************************************************/
#include <stdio.h>
#include <iostream>
#include <arrayfire.h>
using namespace std;
using namespace af;
void simulateHestonModel(af::array &xres, af::array &vres,
float T, unsigned int N, unsigned int R, float mu, float kappa,
float vBar, float sigmaV, float rho, float x0, float v0)
{
float deltaT = T / (float)(N - 1);
af::array x[] = {af::constant(x0, R), af::constant(0, R)};
af::array v[] = {af::constant(v0, R), af::constant(0, R)};
float sqrtDeltaT = sqrt(deltaT);
float sqrtOneMinusRhoSquare = sqrt(1 - rho*rho);
float mArray[] = {rho, sqrtOneMinusRhoSquare};
af::array m(2, 1, mArray);
unsigned int tPrevious = 0, tCurrent = 0;
af::array zeroConstant = constant(0, R);
for (unsigned int t = 1; t < N; t++) {
tPrevious = (t+1) % 2;
tCurrent = t % 2;
af::array dBt = randn(R, 2) * sqrtDeltaT;
af::array sqrtVLag = af::sqrt(v[tPrevious]);
x[tCurrent]= x[tPrevious] + (mu - 0.5 * v[tPrevious]) * deltaT + (sqrtVLag * dBt(span, 0));
af::array vTmp = v[tPrevious] + kappa * (vBar - v[tPrevious]) * deltaT + sigmaV * (sqrtVLag * matmul(dBt, m));
v[tCurrent] = max(vTmp, zeroConstant);
}
xres = x[tCurrent];
vres = v[tCurrent];
}
int main()
{
float T = 1;
unsigned int nT = 10 * T;
unsigned int R_first_run = 1000;
unsigned int R = 20000000;
float x0 = 0; // initial log stock price
float v0 = pow(0.087, 2); // initial volatility
float r = log(1.0319); // risk-free rate
float rho = -0.82; // instantaneous correlation between Brownian motions
float sigmaV = 0.14; // variance of volatility
float kappa = 3.46; // mean reversion speed
float vBar = 0.008; // mean variance
float k = log(0.95); // strike price
// Price European call option
try {
// first run
simulateHestonModel(x, v, T, nT, R_first_run, r, kappa, vBar, sigmaV, rho, x0, v0);
af::sync(); // Ensure the first run is finished
timer::start();
simulateHestonModel(x, v, T, nT, R, r, kappa, vBar, sigmaV, rho, x0, v0);
cout << "Time in simulation: " << timer::stop() << endl;
af::array K = exp(constant(k, x.dims()));
af::array zeroConstant = constant(0, x.dims());
af::array C_CPU = exp(-r * T) * mean(af::max(af::exp(x) - K, zeroConstant));
af_print(C_CPU);
return 0;
} catch (af::exception& e) {
fprintf(stderr, "%s\n", e.what());
return 1;
}
}