Technology

การจำลองการเคลื่อนที่แบบบราวเนียนด้วย Python

2025-05-30 06:53:26

ในบทความนี้เราจะสำรวจการจำลองการเคลื่อนที่แบบบราวเนียน ซึ่งเป็นหนึ่งในแนวคิดพื้นฐานที่สุดในการกำหนดราคาอนุพันธ์ การเคลื่อนที่แบบบราวเนียนเป็นแบบจำลองทางคณิตศาสตร์ที่ใช้จำลองพฤติกรรมของราคาสินทรัพย์เพื่อวัตถุประสงค์ในการกำหนดราคาในสัญญาออปชั่น

วิธีการทั่วไปในการกำหนดราคาตัวเลือกเหล่านี้บนสินทรัพย์คือการจำลองเส้นทางสินทรัพย์แบบสุ่มจำนวนมากตลอดอายุของตัวเลือก กำหนดราคาของตัวเลือกในแต่ละสถานการณ์เหล่านี้ เฉลี่ยราคานี้ แล้วหักส่วนลดเพื่อผลิตราคาสุดท้าย นี่คือตัวอย่างของวิธีมอนติคาร์โล

บทความนี้จะมุ่งเน้นไปที่การจำลองเส้นทางการเคลื่อนที่แบบบราวน์ด้วยภาษาโปรแกรม Python โดยเฉพาะ มันจะเริ่มต้นด้วยการดูการเคลื่อนที่แบบบราวน์มาตรฐาน จากนั้นจะเพิ่มความซับซ้อนให้กับพลศาสตร์ของการเคลื่อนที่แบบบราวน์ เช่นเคยกับบทความของ QuantStart ทุกบทความจะมีโค้ดเต็มเพื่อทำซ้ำผลลัพธ์ที่แน่นอนด้านล่างนี้

บทความนี้ได้รับอิทธิพลอย่างมากจากหนังสือที่แนะนำอย่างยิ่งเรื่อง Monte Carlo Methods for Financial Engineering โดย Paul Glasserman [Glasserman, 2003] หนังสือเล่มนี้ครอบคลุมรายละเอียดทางทฤษฎีที่จำเป็นซึ่งนำไปสู่สูตรการแบ่งส่วนที่นำเสนอ (โดยไม่มีการอธิบายการอนุมาน) ด้านล่างนี้

บทความจะเริ่มต้นด้วยการจำลองการเคลื่อนที่แบบบราวน์ "มาตรฐาน" ซึ่งเกี่ยวข้องกับเส้นทางการเคลื่อนที่แบบบราวน์ที่มีค่าเฉลี่ยเป็นศูนย์และความแปรปรวนเป็นหนึ่ง จากนั้นจะพูดคุยเกี่ยวกับวิธีการรวมค่าเฉลี่ยคงที่ที่ไม่เป็นศูนย์และความแปรปรวนคงที่ที่ไม่เป็นหนึ่งสำหรับการจำลองเส้นทางการเคลื่อนที่แบบบราวเนียน ในบทความต่อไปเราจะพูดคุยเกี่ยวกับพลศาสตร์ของค่าเฉลี่ยและความแปรปรวนที่เปลี่ยนแปลงตามเวลา โดยเริ่มต้นด้วยค่าเฉลี่ยและความแปรปรวนที่คงที่เป็นช่วงๆ และต่อมาเป็นการจำลองด้วยค่าเฉลี่ยและความแปรปรวนที่เปลี่ยนแปลงตามเวลาอย่างเต็มที่

การเคลื่อนที่แบบบราวเนียนมาตรฐาน

สมการต่อไปนี้ถูกนำมาจากสมการ 3.2 ใน [Glasserman, 2003] มันให้ความสัมพันธ์การวนซ้ำสำหรับเส้นทางการเคลื่อนที่แบบบราวน์ (เดี่ยว) ข้ามจุด n จุด (ไม่จำเป็นต้องกระจายเท่ากันในเวลา) ความสัมพันธ์นี้กล่าวโดยสรุปว่าค่าของเส้นทางการเคลื่อนที่แบบบราวเนียนใหม่ที่เวลา $t_{i + 1}$ เท่ากับค่าที่เวลา $t_{i}$

(เช่น ค่าหนึ่งก่อนหน้า) บวกกับรากที่สองของความแตกต่างของเวลาระหว่างแต่ละค่า คูณด้วยตัวอย่างที่แจกแจงตามปกติที่เวลา $t_{i + 1}$

$\begin{array}{r} W (t_{i + 1}) = W (t_{i}) + \sqrt{t_{i + 1} - t_{i}} Z_{i + 1} \end{array}$

อย่างไรก็ตาม เราจะใช้วิธีการเวกเตอร์ของไลบรารี Python NumPy เพื่อสร้างเส้นทาง k พร้อมกัน เพื่อให้บรรลุเป้าหมายนี้ เราจะสุ่มตัวอย่าง $k \times n$ ตัวจากการแจกแจงปกติมาตรฐาน จากนั้นใช้ตัวอย่างเหล่านี้ใน k

การจำลองเส้นทางแบบเวกเตอร์ผ่านการใช้ความสัมพันธ์การวนซ้ำข้างต้น

งานแรกคือการนำเข้าห้องสมุดเชิงปริมาณและการสร้างภาพที่จำเป็นสำหรับการวิเคราะห์ในภายหลัง ซึ่งรวมถึง Matplotlib, NumPy, Pandas และ Seaborn:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

เรายังต้องตั้งค่า seed แบบสุ่มเพื่อให้แน่ใจว่าค่าของคุณตรงกับของเราอย่างแม่นยำ สำหรับสิ่งนี้ เราจำเป็นต้องใช้ตัวสร้างหมายเลขสุ่ม (pseudo-)Random Number Generator (RNG) ซึ่งในกรณีนี้คือค่าเริ่มต้นที่ให้โดย NumPy จากนั้นจำเป็นต้องระบุค่า seed ที่เป็นจำนวนเต็มเพื่อให้ผลลัพธ์ทั้งหมดด้านล่างนี้สามารถทำซ้ำได้:

rng = np.random.default_rng(42)

ขั้นตอนถัดไปคือการเลือกจำนวนเส้นทางและจำนวนช่วงเวลาในแต่ละเส้นทาง เราได้กำหนด k = 50 เส้นทาง โดยแต่ละเส้นทางมี n = 1000 ขั้นตอนเวลา:

paths = 50
points = 1000

เนื่องจากนี่คือการเคลื่อนที่แบบบราวเนียนมาตรฐาน การลอยตัว (ค่าเฉลี่ย) ของกระบวนการนี้เป็นศูนย์ ในขณะที่ความผันผวน (ส่วนเบี่ยงเบนมาตรฐาน) เป็นหนึ่ง:

mu, sigma = 0.0, 1.0

ขั้นตอนถัดไปคือการสุ่มตัวอย่างที่จำเป็นทั้งหมดจากการแจกแจงปกติมาตรฐาน (Gaussian) และเก็บไว้ในเมทริกซ์ที่มีแถวเป็นเส้นทางและคอลัมน์เป็นจุด

( $k \times n$ ค่า). สำหรับสิ่งนี้ เราใช้ RNG เดียวกันกับที่กล่าวถึงข้างต้น เนื่องจากมันถูกตั้งค่าเพื่อให้สามารถทำซ้ำได้:

Z = rng.normal(mu, sigma, (paths, points))

ขั้นตอนถัดไปคือการกำหนดช่วงเวลาที่การจำลองจะเกิดขึ้น เราจะใช้ช่วงเวลา [0.0, 1.0] เราจะจำกัดตัวเองด้วย สำหรับวัตถุประสงค์ของบทแนะนำนี้ ให้ใช้ขนาดก้าวที่เท่ากัน ซึ่งแสดงด้วย $Δ t = t_{i + 1} - t_{i}$ สำหรับสิ่งนี้ เราเพียงแค่หารความยาวของช่วงด้วยจำนวนจุดเวลาแล้วลบหนึ่ง:

interval = [0.0, 1.0]
dt = (interval[1] - interval[0]) / (points - 1)

ในการสร้างกราฟใน Matplotlib จำเป็นต้องสร้างแกน "linear space" ต่อไปนี้จะใช้ช่วงเวลาและสร้างจุดจำนวนที่มีระยะห่างเท่ากัน (รวมถึงค่าต้นและค่าปลายของช่วงเวลา) ในอาร์เรย์ ซึ่งจะใช้เป็น "แกนเวลา" ของเรา:

t_axis = np.linspace(interval[0], interval[1], points)

โค้ดต่อไปนี้ใช้สูตรการวนซ้ำ (3.2) จาก [Glasserman, 2003] เราสร้างแมทริกซ์ขนาด $k \times n$ ชื่อ W เพื่อเก็บผลลัพธ์ของเส้นทางการเคลื่อนที่แบบบราวเนียนที่เต็มไปด้วยศูนย์ จากนั้นเราจะวนลูปตามจำนวนจุดในแกนเวลา สำหรับแต่ละจุดเหล่านี้ เราใช้การเขียนแบบเวกเตอร์ของ NumPy (W[:, real_idx]) เพื่อกำหนดค่าของเส้นทางถัดไป สำหรับแต่ละเส้นทาง ให้เท่ากับค่าของเส้นทางก่อนหน้านี้บวกกับรากที่สองของช่วงเวลา $Δ t$ คูณด้วยตัวอย่างปกติมาตรฐานที่เหมาะสมภายในเมทริกซ์ Z:

W = np.zeros((paths, points))
for idx in range(points - 1):
    real_idx = idx + 1
    W[:, real_idx] = W[:, real_idx - 1] + np.sqrt(dt) * Z[:, idx]

เราสามารถใช้ Matplotlib เพื่อวาดเส้นทางทั้งหมดเหล่านี้บนรูปเดียวกันได้ สำหรับสิ่งนี้ เราใช้วิธี subplots ซึ่งช่วยให้เราควบคุมการกำหนดค่าของรูปภาพและแกนได้อย่างละเอียดมากขึ้น เราระบุว่าเราต้องการแถวและคอลัมน์เดียว จากนั้นระบุขนาดของรูปในหน่วยนิ้ว (ที่ความละเอียดเริ่มต้น 100 จุดต่อนิ้ว (dpi)) จากนั้นเราจะวนลูปตามจำนวนเส้นทางและวาดแกนเวลาเทียบกับแถวที่เหมาะสมในเมทริกซ์ W โดยใช้การเขียนแบบเวกเตอร์ของ NumPy (W[path, :]) สุดท้าย เราวาดกราฟ:

fig, ax = plt.subplots(1, 1, figsize=(12, 8))
for path in range(paths):
    ax.plot(t_axis, W[path, :])
ax.set_title("Standard Brownian Motion sample paths")
ax.set_xlabel("Time")
ax.set_ylabel("Asset Value")
plt.show()

สิ่งนี้ผลิตภาพต่อไปนี้:

ค่อนข้างชัดเจนว่าพื้นที่ส่วนใหญ่ของเส้นทางกำลังรวมตัวกันรอบศูนย์ใกล้ t = 1.0 (ตอนสิ้นสุดของการจำลอง) อยู่ภายในค่าระหว่าง -1 และ 1 ในขณะที่มีตัวอย่างบางอย่างที่มีค่าปลายทางที่สุดขั้วมากกว่า นี่คือพฤติกรรมที่คาดหวังสำหรับการเคลื่อนที่แบบบราวน์มาตรฐาน

อย่างไรก็ตาม สิ่งนี้สามารถทำให้มีความเป็นเชิงปริมาณมากขึ้นโดยการประมาณการแจกแจงของค่าพาธสุดท้ายสำหรับแต่ละ k = 50 เส้นทางตัวอย่างที่ผลิตขึ้น เพื่อให้บรรลุเป้าหมายนี้ เราสามารถใช้ Kernel Density Estimate (KDE) ซึ่งเป็นวิธีการประมาณการการแจกแจงความน่าจะเป็น (ต่อเนื่อง) จากข้อมูลตัวอย่างที่มีขีดจำกัด ห้องสมุดการสร้างภาพทางสถิติ Seaborn จะถูกใช้สำหรับสิ่งนี้

ขั้นตอนแรกคือการสร้าง Pandas DataFrame (ที่มีคอลัมน์เดียว) ของค่าปลายทางของเส้นทางตัวอย่าง:

final_values = pd.DataFrame({'final_values': W[:, -1]})

จากนั้นเราสามารถใช้ Seaborn เพื่อประมาณและสร้างกราฟการกระจาย KDE ของค่าที่ได้สุดท้าย:

fig, ax = plt.subplots(1, 1, figsize=(12, 8))
sns.kdeplot(data=final_values, x='final_values', fill=True, ax=ax)
ax.set_title("Kernel Density Estimate of asset path final value distribution")
ax.set_ylim(0.0, 0.325)
ax.set_xlabel('Final Values of Asset Paths')
plt.show()

สิ่งนี้ผลิตภาพต่อไปนี้:

ดังที่เห็น นี่ดูเหมือนการแจกแจงปกติมาตรฐาน (เช่น ค่าเฉลี่ยเป็นศูนย์ ส่วนเบี่ยงเบนมาตรฐานเป็นหนึ่ง) อย่างไรก็ตาม เราสามารถคำนวณค่าเฉลี่ยและส่วนเบี่ยงเบนมาตรฐานของตัวอย่างจากการแจกแจงนี้ได้เช่นกัน:

print(final_values.mean(), final_values.std())

นี่ผลิตผลลัพธ์ดังต่อไปนี้:

(final_values   -0.011699
 dtype: float64,
 final_values    1.251382
 dtype: float64)

ค่าเฉลี่ย -0.01 ใกล้เคียงกับศูนย์และส่วนเบี่ยงเบนมาตรฐาน 1.25 ค่อนข้างใกล้เคียงกับหนึ่ง หากคำนวณเส้นทางเพิ่มเติม ค่าต่างๆ เหล่านี้จะมีแนวโน้มไปในทิศทางของการแจกแจงปกติมาตรฐาน เนื่องจากกฎของจำนวนมาก

ในส่วนถัดไปเราจะอนุญาตให้มีค่าเฉลี่ยที่ไม่เป็นศูนย์และส่วนเบี่ยงเบนมาตรฐานที่ไม่เป็นหนึ่งสำหรับเส้นทางการสุ่มตัวอย่าง

การล่องลอยและความผันผวนคงที่ การเคลื่อนที่แบบบราวเนียน

ความสัมพันธ์การวนซ้ำถัดไปถูกนำมาจากสมการ 3.3 จาก [Glasserman, 2003] มันเป็นการปรับเปลี่ยนจากสมการข้างต้นเพื่อรวมเทอมเพิ่มเติม $μ (t_{i + 1} - t_{i})$ พร้อมกับการคูณเทอมสุดท้ายด้วย $σ$ นี่คือสิ่งที่ทำให้ค่าเฉลี่ยไม่เป็นศูนย์ รวมถึงส่วนเบี่ยงเบนมาตรฐานที่ไม่เป็นหนึ่ง:

$\begin{array}{r} X (t_{i + 1}) = X (t_{i}) + μ (t_{i + 1} - t_{i}) + σ \sqrt{t_{i + 1} - t_{i}} Z_{i + 1} \end{array}$

ส่วนใหญ่ของโค้ดและตัวแปรที่คำนวณจากส่วนก่อนหน้านี้สามารถนำกลับมาใช้ใหม่ได้ทั้งหมด รวมถึงเมทริกซ์การสุ่มปกติแบบมาตรฐานที่มีการตั้งค่าแบบสุ่ม Z จำเป็นต้องตั้งค่าลอยตัวใหม่ mu_c และความผันผวน sigma_c ซึ่งเราได้ตั้งไว้ที่ 5 และ 2 ตามลำดับ:

mu_c, sigma_c = 5.0, 2.0

โค้ดต่อไปนี้ใช้ในการคำนวณสมการ 3.3 จาก [Glasserman, 2003] สังเกตการรวม mu_c * dt รวมถึงการรวม sigma_c * ... ในส่วนหลัง ยังสามารถใช้ $Δ t = t_{i + 1} - t_{i}$ (dt) ได้เนื่องจากเราใช้ช่วงเวลาคงที่เท่านั้น:

X = np.zeros((paths, points))
for idx in range(points - 1):
    real_idx = idx + 1
    X[:, real_idx] = X[:, real_idx - 1] + mu_c * dt + sigma_c * np.sqrt(dt) * Z[:, idx]

เช่นเดียวกับการเคลื่อนที่แบบบราวน์มาตรฐาน เราสามารถสร้างกราฟเส้นทางตัวอย่างที่ผลิตขึ้นได้:

fig, ax = plt.subplots(1, 1, figsize=(12, 8))
for path in range(paths):
    ax.plot(t_axis, X[path, :])
ax.set_title("Constant mean and standard deviation Brownian Motion sample paths")
ax.set_xlabel("Time")
ax.set_ylabel("Asset Value")
plt.show()

สิ่งนี้ผลิตภาพต่อไปนี้:

สังเกตการเคลื่อนที่เชิงบวกไปยังค่าเฉลี่ยที่ 5 โดยมีการกระจายที่เพิ่มขึ้นเมื่อเปรียบเทียบกับการเคลื่อนที่แบบบราวเนียนมาตรฐาน ยังสามารถเปลี่ยนแปลงค่าของ mu_c และ sigma_c เพื่อสร้างพลศาสตร์ของเส้นทางที่แตกต่างกันได้

ขั้นตอนถัดไป

ในบทความถัดไปในซีรีส์นี้ เราจะผ่อนคลายสมมติฐานของค่าเฉลี่ยและส่วนเบี่ยงเบนมาตรฐานคงที่ โดยอนุญาตให้มีค่าเฉลี่ยและส่วนเบี่ยงเบนมาตรฐานที่เปลี่ยนแปลงตามเวลาเพื่อแนะนำพลศาสตร์ของสินทรัพย์ที่ซับซ้อนมากขึ้น ในบทความต่อไปเราจะพิจารณาการเคลื่อนที่แบบบราวเนียนหลายมิติ "บราวเนียนบริดจ์" และการเคลื่อนที่แบบบราวเนียนเชิงเรขาคณิต

โค้ดทั้งหมด

# Imports
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

# Seed the random number generator
rng = np.random.default_rng(42)

# Determine the number of paths and points per path
points = 1000
paths = 50

# Create the initial set of random normal draws
mu, sigma = 0.0, 1.0
Z = rng.normal(mu, sigma, (paths, points))

# Define the time step size and t-axis
interval = [0.0, 1.0]
dt = (interval[1] - interval[0]) / (points - 1)
t_axis = np.linspace(interval[0], interval[1], points)

# Use Equation 3.2 from [Glasserman, 2003] to sample 50 standard brownian motion paths
W = np.zeros((paths, points))
for idx in range(points - 1):
    real_idx = idx + 1
    W[:, real_idx] = W[:, real_idx - 1] + np.sqrt(dt) * Z[:, idx]

# Plot these paths
fig, ax = plt.subplots(1, 1, figsize=(12, 8))
for path in range(paths):
    ax.plot(t_axis, W[path, :])
ax.set_title("Standard Brownian Motion sample paths")
ax.set_xlabel("Time")
ax.set_ylabel("Asset Value")
plt.show()

# Obtain the set of final path values
final_values = pd.DataFrame({'final_values': W[:, -1]})

# Estimate and plot the distribution of these final values with Seaborn
fig, ax = plt.subplots(1, 1, figsize=(12, 8))
sns.kdeplot(data=final_values, x='final_values', fill=True, ax=ax)
ax.set_title("Kernel Density Estimate of asset path final value distribution")
ax.set_ylim(0.0, 0.325)
ax.set_xlabel('Final Values of Asset Paths')
plt.show()

# Output the mean and stdev of these final values
print(final_values.mean(), final_values.std())

# Create a non-zero mean and non-unit standard deviation
mu_c, sigma_c = 5.0, 2.0

# Use Equation 3.3 from [Glasserman, 2003] to sample 50 brownian motion paths
X = np.zeros((paths, points))
for idx in range(points - 1):
    real_idx = idx + 1
    X[:, real_idx] = X[:, real_idx - 1] + mu_c * dt + sigma_c * np.sqrt(dt) * Z[:, idx]

# Plot these paths
fig, ax = plt.subplots(1, 1, figsize=(12, 8))
for path in range(paths):
    ax.plot(t_axis, X[path, :])
ax.set_title("Constant mean and standard deviation Brownian Motion sample paths")
ax.set_xlabel("Time")
ax.set_ylabel("Asset Value")
plt.show()

อ้างอิง :Brownian Motion Simulation with Python

จาก https://www.quantstart.com/articles/brownian-motion-simulation-with-python/