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 <arrayfire.h>
#include <stdio.h>
#include <iostream>
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;
}
}
af::pow
AFAPI array pow(const array &lhs, const array &rhs)
C++ Interface for power.
af::matmul
AFAPI array matmul(const array &lhs, const array &rhs, const matProp optLhs=AF_MAT_NONE, const matProp optRhs=AF_MAT_NONE)
Matrix multiply of two arrays.
af::exp
AFAPI array exp(const array &in)
C++ Interface for exponential of an array.
af::constant
array constant(T val, const dim4 &dims, const dtype ty=(af_dtype) dtype_traits< T >::ctype)
af::log
AFAPI array log(const array &in)
C++ Interface for natural logarithm.
af::mean
AFAPI array mean(const array &in, const dim_t dim=-1)
C++ Interface for mean.
af::array
A multi dimensional data container.
Definition: array.h:35
af
Definition: algorithm.h:15
af_print
#define af_print(...)
Definition: util.h:148
af::max
AFAPI array max(const array &in, const int dim=-1)
C++ Interface for maximum values in an array.
af::randn
AFAPI array randn(const dim4 &dims, const dtype ty, randomEngine &r)
af::sqrt
AFAPI array sqrt(const array &in)
C++ Interface for square root of input.
af::exception
An ArrayFire exception class.
Definition: exception.h:29
af::span
AFAPI seq span
A special value representing the entire axis of an af::array.
af::array::dims
dim4 dims() const
Get dimensions of the array.
arrayfire.h
af::sync
AFAPI void sync(const int device=-1)
Blocks until the device is finished processing.
af::exception::what
virtual const char * what() const
Returns an error message for the exception in a string format.
Definition: exception.h:60