模拟几何布朗运动
Simulating the Geometric Brownian Motion
背景资料:
考虑伪代码:
问题:
我似乎有问题的是算法 1 伪代码的 S 上标 j 部分,我的教授说由于 N = 10,000 和 n = 10,我正在计算 10 个 S。这是我遇到问题的代码部分:
double Z[n][N];
for(int j = 0; j < N; j++){
for(int i = 0; i < n/2; i++){
Z[2*i][j] = Box_MullerX(n,N,shifted_hs[i],shifted_hs[i+1]);
Z[2*i+1][j] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
}
}
/*for(int j = 0; j < N; j++){
for(int i = 0; i < n; i++){
cout << Z[j][i] << " ";
}
cout << endl;
}*/
double S[N][n];
S[0][0] = S0;
double t[n+1];
for(int i = 0; i <= n; i++){
t[i] = (double)i*T/(n-1);
}
for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
S[j][i] = S[j-1][i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*Z[j][i]);
}
}
具体来说,我只需要计算代码的 S^j 部分的帮助(假设其他一切都正确)。
为了完整起见,这是我的全部代码:
#include <iostream>
#include <cmath>
#include <math.h>
#include <fstream>
#include <random>
using namespace std;
double phi(long double x);
double Black_Scholes(double T, double K, double S0, double r, double sigma);
double Halton_Seq(int index, double base);
double Box_MullerX(int n, int N,double u_1, double u_2);
double Box_MullerY(int n, int N,double u_1, double u_2);
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0);
int main(){
int n = 10;
int N = 10000;
// Calculate BS
double T = 1.0;
double K = 50.0;
double r = 0.1;
double sigma = 0.3;
double S0 = 50.0;
double price = Black_Scholes(T,K,S0,r,sigma);
//cout << price << endl;
// Random-shift Halton Sequence
Random_Shift_Halton(n,N,r,sigma,T, S0);
}
double phi(double x){
return 0.5 * erfc(-x * M_SQRT1_2);
}
double Black_Scholes(double T, double K, double S0, double r, double sigma){
double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double d_2 = (log(S0/K) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double C = S0*phi(d_1) - phi(d_2)*K*exp(-r*T);
return C;
}
double Halton_Seq(int index, double base){
double f = 1.0, r = 0.0;
while(index > 0){
f = f/base;
r = r + f*(fmod(index,base));
index = index/base;
}
return r;
}
double Box_MullerX(int n, int N,double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2*M_PI*u_2;
double X = sqrt(R)*cos(theta);
return X;
}
double Box_MullerY(int n, int N,double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2*M_PI*u_2;
double Y = sqrt(R)*sin(theta);
return Y;
}
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0){
// Generate the Halton_Sequence
double hs[N + 1];
for(int i = 0; i <= N; i++){
hs[i] = Halton_Seq(i,2.0);
}
// Generate two random numbers from U(0,1)
double u[N+1];
random_device rd;
mt19937 e2(rd());
uniform_real_distribution<double> dist(0, 1);
for(int i = 0; i <= N; i++){
u[i] = dist(e2);
}
// Add the vector u to each vector hs and take fraction part of the sum
double shifted_hs[N+1];
for(int i = 0; i <= N; i++){
shifted_hs[i] = hs[i] + u[i];
if(shifted_hs[i] > 1){
shifted_hs[i] = shifted_hs[i] - 1;
}
}
// Use Geometric Brownian Motion
/*double Z[N];
for(int i = 0; i < N; i++){
if(i % 2 == 0){
Z[i] = Box_MullerX(n,N,shifted_hs[i+1],shifted_hs[i+2]);
}else{
Z[i] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
}
}*/
double Z[n][N];
for(int j = 0; j < N; j++){
for(int i = 0; i < n/2; i++){
Z[2*i][j] = Box_MullerX(n,N,shifted_hs[i],shifted_hs[i+1]);
Z[2*i+1][j] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
}
}
/*for(int j = 0; j < N; j++){
for(int i = 0; i < n; i++){
cout << Z[j][i] << " ";
}
cout << endl;
}*/
double S[N][n];
S[0][0] = S0;
double t[n+1];
for(int i = 0; i <= n; i++){
t[i] = (double)i*T/(n-1);
}
for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
S[j][i] = S[j-1][i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*Z[j][i]);
}
}
for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
cout << S[j][i] << endl;
}
}
// Use random-shift halton sequence to obtain 40 independent estimates for the price of the option
}
请说明您的期望和得到的是什么。这有助于理解可能的变体。
而且我注意到您在代码中混合了整数和双精度值。这不是好的做法。
另外 - Mark Setchell 对 - 检查你的 sigma 值。我想它应该是平方的。并尝试正确设置 S
的初始值。
我尝试按照我理解的方式修复您代码中的错误。由于编译器版本较旧,代码略有改动。
#include <iostream>
#include <cmath>
#include <math.h>
#include <random>
using namespace std;
const double M_PI = 3.1415926535897932384626433832795;
const double M_SQRT1_2 = 0.70710678118654752440084436210485;
double phi(long double x);
double Black_Scholes(double T, double K, double S0, double r, double sigma);
double Halton_Seq(int index, double base);
double Box_MullerX(double u_1, double u_2);
double Box_MullerY(double u_1, double u_2);
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0);
double erfc(double v) {
double y = 1.0 / ( 1.0 + 0.3275911 * (v >= 0 ? v : -v));
y = 1 - (((((
+ 1.061405429 * y
- 1.453152027) * y
+ 1.421413741) * y
- 0.284496736) * y
+ 0.254829592) * y)
* exp (-v * v);
return (v >= 0 ? y : -y);
}
double phi(double x){
return 0.5 * erfc(-x * M_SQRT1_2);
}
double Black_Scholes(double T, double K, double S0, double r, double sigma){
double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double d_2 = (log(S0/K) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double C = S0*phi(d_1) - phi(d_2)*K*exp(-r*T);
return C;
}
double Halton_Seq(int index, double base){
double f = 1.0, r = 0.0;
double indice = (double)index;
while(indice >= base){
f = f/base;
r = r + f*(fmod(indice,base));
indice = ceil(indice/base);
}
return r;
}
double Box_MullerX(double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2.0*M_PI*u_2;
double X = sqrt(R)*cos(theta);
return X;
}
double Box_MullerY(double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2.0*M_PI*u_2;
double Y = sqrt(R)*sin(theta);
return Y;
}
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0){
// Generate the Halton_Sequence
double* hs = new double[N+1];
for(int i = 0; i <= N; i++){
hs[i] = Halton_Seq(i,2.0);
}
// Generate two random numbers from U(0,1)
double* u = new double[N+1];
random_device rd;
mt19937 e2(rd());
uniform_real_distribution<double> dist(0.0, 1.0);
for(int i = 0; i <= N; i++){
u[i] = dist(e2);
}
// Add the vector u to each vector hs and take fraction part of the sum
double* shifted_hs = new double[N+1];
for(int i = 0; i <= N; i++){
shifted_hs[i] = hs[i] + u[i];
if(shifted_hs[i] > 1){
shifted_hs[i] = shifted_hs[i] - 1.0;
}
}
// Use Geometric Brownian Motion
double** Z = new double*[N];
for(int i = 0; i < N; ++i)
Z[i] = new double[n];
for(int j = 0; j < N; j++){
for(int i = 0; i < n/2; i++){
Z[j][2*i] = Box_MullerX(shifted_hs[j],shifted_hs[j+1]);
Z[j][2*i+1] = Box_MullerY(shifted_hs[j],shifted_hs[j+1]);
}
}
double** S = new double*[N+1];
for(int i = 0; i <= N; ++i)
S[i] = new double[n+1];
for (int i = 0; i <= N; i++) S[i][0] = S0;
for (int j = 0; j <= n; j++) S[0][j] = S0;
double* t = new double[n+1];
for(int i = 0; i <= n; i++){
t[i] = (double)i*T/(double)(n-1);
}
for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
S[j][i] = S[j-1][i-1]*exp((mu - sigma/2.0)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*Z[j-1][i-1]);
}
}
// Use random-shift halton sequence to obtain 40 independent estimates for the price of the option
//....
delete(hs);
delete(u);
delete(shifted_hs);
for(int i = 0; i < N; ++i) {
delete(Z[i]);
}
for(int i = 0; i <= N; ++i) delete(S[i]);
delete(S);
delete(Z);
delete(t);
return 0.0;
}
int main()
{
int n = 10;
int N = 10000;
// Calculate BS
double T = 1.0;
double K = 50.0;
double r = 0.1;
double sigma = 0.3;
sigma *= sigma;
double S0 = 50.0;
double price = Black_Scholes(T,K,S0,r,sigma);
//cout << price << endl;
// Random-shift Halton Sequence
Random_Shift_Halton(n,N,r,sigma,T, S0);
return 0;
}
这就是你要找的吗?
背景资料:
考虑伪代码:
问题:
我似乎有问题的是算法 1 伪代码的 S 上标 j 部分,我的教授说由于 N = 10,000 和 n = 10,我正在计算 10 个 S。这是我遇到问题的代码部分:
double Z[n][N];
for(int j = 0; j < N; j++){
for(int i = 0; i < n/2; i++){
Z[2*i][j] = Box_MullerX(n,N,shifted_hs[i],shifted_hs[i+1]);
Z[2*i+1][j] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
}
}
/*for(int j = 0; j < N; j++){
for(int i = 0; i < n; i++){
cout << Z[j][i] << " ";
}
cout << endl;
}*/
double S[N][n];
S[0][0] = S0;
double t[n+1];
for(int i = 0; i <= n; i++){
t[i] = (double)i*T/(n-1);
}
for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
S[j][i] = S[j-1][i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*Z[j][i]);
}
}
具体来说,我只需要计算代码的 S^j 部分的帮助(假设其他一切都正确)。
为了完整起见,这是我的全部代码:
#include <iostream>
#include <cmath>
#include <math.h>
#include <fstream>
#include <random>
using namespace std;
double phi(long double x);
double Black_Scholes(double T, double K, double S0, double r, double sigma);
double Halton_Seq(int index, double base);
double Box_MullerX(int n, int N,double u_1, double u_2);
double Box_MullerY(int n, int N,double u_1, double u_2);
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0);
int main(){
int n = 10;
int N = 10000;
// Calculate BS
double T = 1.0;
double K = 50.0;
double r = 0.1;
double sigma = 0.3;
double S0 = 50.0;
double price = Black_Scholes(T,K,S0,r,sigma);
//cout << price << endl;
// Random-shift Halton Sequence
Random_Shift_Halton(n,N,r,sigma,T, S0);
}
double phi(double x){
return 0.5 * erfc(-x * M_SQRT1_2);
}
double Black_Scholes(double T, double K, double S0, double r, double sigma){
double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double d_2 = (log(S0/K) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double C = S0*phi(d_1) - phi(d_2)*K*exp(-r*T);
return C;
}
double Halton_Seq(int index, double base){
double f = 1.0, r = 0.0;
while(index > 0){
f = f/base;
r = r + f*(fmod(index,base));
index = index/base;
}
return r;
}
double Box_MullerX(int n, int N,double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2*M_PI*u_2;
double X = sqrt(R)*cos(theta);
return X;
}
double Box_MullerY(int n, int N,double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2*M_PI*u_2;
double Y = sqrt(R)*sin(theta);
return Y;
}
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0){
// Generate the Halton_Sequence
double hs[N + 1];
for(int i = 0; i <= N; i++){
hs[i] = Halton_Seq(i,2.0);
}
// Generate two random numbers from U(0,1)
double u[N+1];
random_device rd;
mt19937 e2(rd());
uniform_real_distribution<double> dist(0, 1);
for(int i = 0; i <= N; i++){
u[i] = dist(e2);
}
// Add the vector u to each vector hs and take fraction part of the sum
double shifted_hs[N+1];
for(int i = 0; i <= N; i++){
shifted_hs[i] = hs[i] + u[i];
if(shifted_hs[i] > 1){
shifted_hs[i] = shifted_hs[i] - 1;
}
}
// Use Geometric Brownian Motion
/*double Z[N];
for(int i = 0; i < N; i++){
if(i % 2 == 0){
Z[i] = Box_MullerX(n,N,shifted_hs[i+1],shifted_hs[i+2]);
}else{
Z[i] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
}
}*/
double Z[n][N];
for(int j = 0; j < N; j++){
for(int i = 0; i < n/2; i++){
Z[2*i][j] = Box_MullerX(n,N,shifted_hs[i],shifted_hs[i+1]);
Z[2*i+1][j] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
}
}
/*for(int j = 0; j < N; j++){
for(int i = 0; i < n; i++){
cout << Z[j][i] << " ";
}
cout << endl;
}*/
double S[N][n];
S[0][0] = S0;
double t[n+1];
for(int i = 0; i <= n; i++){
t[i] = (double)i*T/(n-1);
}
for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
S[j][i] = S[j-1][i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*Z[j][i]);
}
}
for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
cout << S[j][i] << endl;
}
}
// Use random-shift halton sequence to obtain 40 independent estimates for the price of the option
}
请说明您的期望和得到的是什么。这有助于理解可能的变体。
而且我注意到您在代码中混合了整数和双精度值。这不是好的做法。
另外 - Mark Setchell 对 - 检查你的 sigma 值。我想它应该是平方的。并尝试正确设置 S
的初始值。
我尝试按照我理解的方式修复您代码中的错误。由于编译器版本较旧,代码略有改动。
#include <iostream>
#include <cmath>
#include <math.h>
#include <random>
using namespace std;
const double M_PI = 3.1415926535897932384626433832795;
const double M_SQRT1_2 = 0.70710678118654752440084436210485;
double phi(long double x);
double Black_Scholes(double T, double K, double S0, double r, double sigma);
double Halton_Seq(int index, double base);
double Box_MullerX(double u_1, double u_2);
double Box_MullerY(double u_1, double u_2);
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0);
double erfc(double v) {
double y = 1.0 / ( 1.0 + 0.3275911 * (v >= 0 ? v : -v));
y = 1 - (((((
+ 1.061405429 * y
- 1.453152027) * y
+ 1.421413741) * y
- 0.284496736) * y
+ 0.254829592) * y)
* exp (-v * v);
return (v >= 0 ? y : -y);
}
double phi(double x){
return 0.5 * erfc(-x * M_SQRT1_2);
}
double Black_Scholes(double T, double K, double S0, double r, double sigma){
double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double d_2 = (log(S0/K) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double C = S0*phi(d_1) - phi(d_2)*K*exp(-r*T);
return C;
}
double Halton_Seq(int index, double base){
double f = 1.0, r = 0.0;
double indice = (double)index;
while(indice >= base){
f = f/base;
r = r + f*(fmod(indice,base));
indice = ceil(indice/base);
}
return r;
}
double Box_MullerX(double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2.0*M_PI*u_2;
double X = sqrt(R)*cos(theta);
return X;
}
double Box_MullerY(double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2.0*M_PI*u_2;
double Y = sqrt(R)*sin(theta);
return Y;
}
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0){
// Generate the Halton_Sequence
double* hs = new double[N+1];
for(int i = 0; i <= N; i++){
hs[i] = Halton_Seq(i,2.0);
}
// Generate two random numbers from U(0,1)
double* u = new double[N+1];
random_device rd;
mt19937 e2(rd());
uniform_real_distribution<double> dist(0.0, 1.0);
for(int i = 0; i <= N; i++){
u[i] = dist(e2);
}
// Add the vector u to each vector hs and take fraction part of the sum
double* shifted_hs = new double[N+1];
for(int i = 0; i <= N; i++){
shifted_hs[i] = hs[i] + u[i];
if(shifted_hs[i] > 1){
shifted_hs[i] = shifted_hs[i] - 1.0;
}
}
// Use Geometric Brownian Motion
double** Z = new double*[N];
for(int i = 0; i < N; ++i)
Z[i] = new double[n];
for(int j = 0; j < N; j++){
for(int i = 0; i < n/2; i++){
Z[j][2*i] = Box_MullerX(shifted_hs[j],shifted_hs[j+1]);
Z[j][2*i+1] = Box_MullerY(shifted_hs[j],shifted_hs[j+1]);
}
}
double** S = new double*[N+1];
for(int i = 0; i <= N; ++i)
S[i] = new double[n+1];
for (int i = 0; i <= N; i++) S[i][0] = S0;
for (int j = 0; j <= n; j++) S[0][j] = S0;
double* t = new double[n+1];
for(int i = 0; i <= n; i++){
t[i] = (double)i*T/(double)(n-1);
}
for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
S[j][i] = S[j-1][i-1]*exp((mu - sigma/2.0)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*Z[j-1][i-1]);
}
}
// Use random-shift halton sequence to obtain 40 independent estimates for the price of the option
//....
delete(hs);
delete(u);
delete(shifted_hs);
for(int i = 0; i < N; ++i) {
delete(Z[i]);
}
for(int i = 0; i <= N; ++i) delete(S[i]);
delete(S);
delete(Z);
delete(t);
return 0.0;
}
int main()
{
int n = 10;
int N = 10000;
// Calculate BS
double T = 1.0;
double K = 50.0;
double r = 0.1;
double sigma = 0.3;
sigma *= sigma;
double S0 = 50.0;
double price = Black_Scholes(T,K,S0,r,sigma);
//cout << price << endl;
// Random-shift Halton Sequence
Random_Shift_Halton(n,N,r,sigma,T, S0);
return 0;
}
这就是你要找的吗?