程式

使用隨機移位 Halton 序列獲得歐洲看漲期權價格的 40 個獨立估計

  • March 11, 2017

背景資料：

隨機移位 Halton 序列：考慮維度中的前六個 Halton 向量 $ 2 $ , 使用基 $ 2 $ 和 $ 3 $ :

$$ \begin{bmatrix} 1/2\ 1/3 \end{bmatrix}, \begin{bmatrix} 1/4\ 2/3 \end{bmatrix},\begin{bmatrix} 3/4\ 1/9 \end{bmatrix},\begin{bmatrix} 1/8\ 4/9 \end{bmatrix},\begin{bmatrix} 5/8\ 7/9 \end{bmatrix}, \begin{bmatrix} 3/8\ 2/9 \end{bmatrix} $$ 要獲得該集合的獨立隨機化，請生成兩個隨機數 $ u_1 $ , $ u_2 $ 從 $ U(0,1) $ ， 說， $ u_1 = 0.7 $ 和 $ u_2 = 0.1 $ . 添加向量 $ u = \begin{bmatrix} u_1\ u_2 \end{bmatrix} $ 到上面的每個向量，按分量，並取總和的小數部分。例如，第一個向量變為

$$ \begin{bmatrix} {1/2 + 0.7}\ {1/3 + 0.1} \end{bmatrix}= \begin{bmatrix} 0.2\ 0.4\overline{3} \end{bmatrix} $$ （這裡 $ {x} $ 表示小數部分 $ x $ ）。同樣，添加 $ u $ 到其餘的向量。這是模擬幾何布朗運動模型的虛擬碼： 在此處輸入圖像描述 問題：

考慮一個帶有參數的歐式看漲期權： $ T = 1 $ , $ K = 50 $ , $ r = 0.1 $ , $ \sigma = 0.3 $ ， 和 $ S_0 = 50 $

使用隨機移位 Halton 序列獲得 $ 40 $ 對期權價格的“獨立”估計。對於每個估計，使用 $ N = 10,000 $ 價格路徑。要模擬路徑，您將模擬幾何布朗運動模型 $ \mu = r $ , 並使用 $ 10 $ 時間步長 $ t_0 = 0 $ , $ t_1 = \Delta t $ , $ t_2 = 2\Delta t,\ldots,t_{10} = 10\Delta t = T $ . 使用 Box-Muller 方法生成標準法線數。

我已經能夠成功地使用 Box-Muller 方法來生成我檢查正確的標準法線數。我還模擬了幾何布朗運動模型 $ \mu = r $ 但我不完全確定它是否正確。

我感到困惑的部分是假設幾何布朗運動模型是正確的，我如何應用隨機移位 Halton 序列來獲得 $ 40 $ 對期權價格的“獨立”估計。

我用 C++ 編寫了後者。這是應該很容易閱讀的程式碼：

#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 ZZ[N + 1];
   ZZ[0] = 0.0;
   for(int i = 1; i <= N; i++){
           ZZ[i] = Z[i-1];
   }

   for(int i = 1; i <= N; i++){
       //cout << ZZ[i] << endl;
   }

   double S[N + 1];
   S[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[i] = S[i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*ZZ[i]);
       }
   }

   for(int i = 1; i <= n; i++){
       //cout << S[i] << endl;
   }

   // Use random-shift halton sequence to obtain 40 independent estimates for the price of the option

}

}

如果有人對我從這裡做的事情有任何建議，我將不勝感激。

實際問題

我想這是作業，所以我只會概述步驟。我理解這個問題的方式如下：

  1. 建構一個函式，模擬標的資產的 10,000 個不同樣本路徑，每個路徑具有 10 個等距時間步長。
  2. 對於這些路徑中的每一個，計算終端歐式看漲期權收益。他們的貼現樣本平均值是您對歐洲看漲期權價格的估計。

對 Halton 序列使用不同的隨機偏移量重複步驟 1 和 2 總共 40 次，以獲得 40 個不同的估計值並記錄結果。

其他程式碼

查看您的程式碼，我可以看到很多問題。以下是一些幫助您入門的提示：

  • 您根本不執行 Halton 序列的隨機化。這應該發生在生成標準 Halton 序列和在 Box-Muller 變換中使用它們以生成隨機法線之間。您發布的摘錄非常清楚它是如何完成的。
  • 此程式碼片段
for(int i = 0; i < N; i++){
   if(i % 2 == 0){
       Z[i] = Box_MullerX(n,N,hs[i+1],hs[i+2]);
   }else{
       Z[i] = Box_MullerY(n,N,hs[i],hs[i+1]);
   }
}

只是碰巧沒有越界，因為N它甚至在你的主線中。

  • 如果您想生成N樣本路徑，那麼如果您只想跟踪每條路徑的目前位置（這已經足夠了） ，那麼您就需要S維度。N即S[j]是目前點 $ j $ -th 路徑，您必須將所有條目初始化S為S0.
  • 考慮您的程式碼片段
for(int j = 1; j <= N; j++){
   for(int i = 1; i <= n; i++){
       S[i] = S[i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*ZZ[i]);
   }
}

除了加熱你的 CPU 外，外循環什麼​​都不做。您從不使用循環變數j，而只是執行完全相同的內部循環N時間。換句話說，您只模擬單個路徑而不是N.

  • Box-Muller 變換需要在單位平方上獨立統一。您使用來自同一個 Halton 序列的兩個連續項目。如果您創建所謂的正常變數的散點圖，您會得到類似於下面第一個圖的內容。第二個是使用兩個具有不同鹼基的 Halton 序列產生的。另見 Jaeckel (2002) 第 9.3.1 章中的相應討論。

Halton Box-Muller 出錯了 Halton Box-Muller 做得對

  • 參考

Jaeckel, Peter (2002)金融工程中的蒙地卡羅方法： John Wiley & Sons

引用自：https://quant.stackexchange.com/questions/32942

相關問答

程式

使用原始蒙地卡羅

  • April 3, 2017
檢視問答
程式

美式期權定價的前向蒙特卡羅方法

  • March 23, 2017
檢視問答
程式

兩種資產最大值的看漲期權的蒙地卡羅近似

  • September 9, 2021
檢視問答
程式

使用 sobol 序列的準隨機標準正態數在 Python 中的蒙地卡羅模擬給出了錯誤的值

  • January 30, 2021
檢視問答
程式

如何解釋和定義 GBM 輸出的統計數據

  • January 12, 2021
檢視問答
程式

Python 蒙地卡羅收斂

  • September 22, 2020
檢視問答