1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
/**********************************************************************************************
* 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 {
af::array x;
af::array v;
// 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);
af::sync();
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;
}
}
|
