Introduction to Numerical Methods

เริ่มจากราคาที่เราต้องการประมาณ

ภายใต้สมมติฐาน no-arbitrage ของ Black–Scholes เราใช้กระบวนการราคาภายใต้ risk-neutral measure $\mathbb Q$

ราคา ณ เวลา t ของสัญญาที่จ่าย $g(S_T)$ เมื่อหมดอายุคือ

สำหรับ Call ใช้ $g(S_T)=\max(S_T-K,0)$ ส่วน Put ใช้ $\max(K-S_T,0)$ ต้องแยก payoff ณ วันหมดอายุ ออกจาก ราคาวันนี้ และกำไรหลังหัก premium

เหตุผลที่ drift เป็น r ไม่ใช่ผลตอบแทนคาดหวังจริง μ มาจากการตีราคาที่สอดคล้องกับพอร์ตเลียนแบบและ no-arbitrage ไม่ใช่การทำนายว่าหุ้นจริงจะโตเท่าดอกเบี้ย หรือการสมมติว่าผู้ลงทุนทุกคนไม่กลัวความเสี่ยง ทบทวนที่ พอร์ตเลียนแบบและ risk-neutral pricing

สำหรับตัวอย่างนี้ สูตร Black–Scholes ให้ Call ประมาณ 9.4134 ดอลลาร์ และ Put ประมาณ 6.4580 ดอลลาร์ เราจะใช้เป็น benchmark ของทั้งสองวิธี ไม่ใช่ราคาตลาดที่สังเกตมา

สุ่มราคาปลายทางในครั้งเดียว

จาก Itô’s lemma สำหรับ log S เราได้

นี่คือ exact terminal simulation ภายใต้ GBM ที่พารามิเตอร์คงที่ จึงไม่ต้องแบ่งวันเพื่อหาราคา European Call/Put ที่ payoff ขึ้นกับราคาปลายทางเพียงค่าเดียว คำว่า exact หมายถึงไม่มีความคลาดเคลื่อนจากการแบ่งเวลาในสูตรนี้ การเฉลี่ยจากจำนวนตัวอย่างจำกัดยังมีความคลาดเคลื่อนอยู่

เมื่อแทนค่าตั้งต้น ถ้าช็อกหนึ่งรอบเป็น Z=0 จะได้ $S_T=100e^{0.01}\approx101.0050$ และ discounted Call payoff ประมาณ 0.9753 ดอลลาร์ นี่คือผลจาก หนึ่งช็อก ไม่ใช่ราคา Call ทั้งสัญญา

ขั้นตอนคำนวณมีดังนี้

1. สุ่ม $Z_j$ จาก Standard Normal อย่างอิสระ
2. คำนวณ $S_T^{(j)}$ ด้วยสูตร exponential
3. คำนวณ discounted payoff $Y_j=e^{-rT}g(S_T^{(j)})$
4. ทำซ้ำ N รอบ แล้วเฉลี่ย $\widehat V_N=N^{-1}\sum_{j=1}^N Y_j$

ถ้า payoff ขึ้นกับเส้นทางล่ะ? Asian Option ที่ใช้ค่าเฉลี่ยราคาตามวันสังเกตต้องจำลองราคาทุกวันที่สัญญาระบุ เรายังใช้ exact GBM update ทีละช่วงได้ แต่การสุ่มเพียง $S_T$ ไม่ให้ข้อมูลระหว่างทาง สำหรับ barrier ที่เฝ้าราคาต่อเนื่อง แม้จำลองจุดกริดได้ตรงตาม GBM ก็ยังอาจพลาดการข้าม barrier ระหว่างจุด

หาก SDE ไม่มีวิธีจำลอง exact ที่สะดวก อาจใช้ Euler–Maruyama

สูตร Euler ในระดับราคาอาจให้ราคาติดลบ และมี discretization error ภายใต้เงื่อนไขมาตรฐานที่เหมาะสม weak error ของค่าคาดหมายมีอันดับ $O(\Delta t)$ ส่วน strong error ที่เทียบเส้นทางด้วย Brownian motion เดียวกันโดยทั่วไปเป็น $O(\sqrt{\Delta t})$ อัตราเหล่านี้ต้องตรวจเงื่อนไขของ SDE และ payoff ไม่ใช่กฎที่ใช้ได้กับทุกสัญญา Mike Giles, Lecture 9

เอกสารต้นทางเสนอ $\sum_{i=1}^{12}U_i-6$ เป็นวิธีประมาณ Normal อย่างง่าย ค่านี้มี mean 0 และ variance 1 แต่มีช่วงจำกัด [−6,6] จึงไม่ใช่ Normal จริง ห้องทดลองใช้ Box–Muller แปลง uniform pseudo-random numbers เป็น Normal แทน และระบุ seed เพื่อให้ทำซ้ำได้ วิธี Box–Muller

ราคาเฉลี่ยต้องมาพร้อมความคลาดเคลื่อน

ให้ $s_Y$ เป็น sample standard deviation ของ discounted payoffs ไม่ใช่ของราคาหุ้น

Standard error วัดการแกว่งของค่าประมาณราคาเมื่อสุ่มใหม่ ช่วงความเชื่อมั่นโดยประมาณ 95% จาก Normal approximation คือ

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

ช่วงนี้อธิบายความไม่แน่นอนของ ค่าประมาณราคาในแบบจำลอง ไม่ใช่ช่วงที่ราคาหุ้นจะอยู่ ไม่ใช่ขอบเขตขาดทุนสูงสุด และไม่รวม model error หากสุ่มซ้ำหลายชุด ช่วง 95% ก็ไม่จำเป็นต้องครอบ benchmark ทุกชุด

เมื่อเพิ่ม N เป็น 4 เท่า SE จะลดลงประมาณครึ่งหนึ่งในระยะยาว เช่น ถ้า SE≈0.14 ดอลลาร์ที่ 10,000 รอบ จะต้องใช้ประมาณ 40,000 รอบเพื่อให้เหลือ 0.07 ดอลลาร์ ค่า error ที่เกิดขึ้นจริงของการรันหนึ่งครั้งไม่จำเป็นต้องลดลงทุกครั้งที่เพิ่ม N

แยกความคลาดเคลื่อนอย่างน้อยสามชั้นก่อนตัดสินใจเพิ่มรอบคำนวณ

การสุ่มเพิ่มลดชั้นแรก การลด time step แก้ชั้นที่สอง ทั้งสองอย่างไม่ทำให้ volatility คงที่กลายเป็นสมมติฐานที่เหมาะกับทุกตลาด

ทดลอง: เพิ่มจำนวนรอบแล้วราคาเปลี่ยนอย่างไร

เริ่มจาก Call ที่ค่าเดิม เพิ่ม N โดยคง seed แล้วสังเกตราคาและ SE จากนั้นกดสุ่มชุดใหม่เพื่อดูการเปลี่ยนแปลงข้ามชุด กราฟใช้แกน log ของ N เพื่อให้มองเห็นทั้งหลักร้อยและหลักหมื่น

ห้องทดลองจำลอง $S_T$ แบบ exact และใช้ pseudo-random generator เดิมของซีรีส์กับ Box–Muller เพื่อการศึกษา Notebook ใช้ตัวสุ่มของ Python จึงไม่จำเป็นต้องได้ราคา Monte Carlo เท่ากับเว็บทุกหลัก แม้ใส่ seed เลขเดียวกัน

อีกมุมหนึ่ง: แทนอนุพันธ์ด้วยค่าบนกริด

Finite difference เริ่มจากการเก็บค่า Option บนจุดราคาและเวลา แทนการสุ่มเส้นทางจำนวนมาก ให้

i วิ่งจาก 0 ถึง M และ k วิ่งจาก 0 ถึง L เราใช้ $\tau=T-t$ คือ เวลาที่เหลือถึงวันหมดอายุ เมื่อเริ่มที่ τ=0 เรารู้ payoff อยู่แล้ว จากนั้นเพิ่ม τ จนถึง T เพื่อหาค่าที่วันนี้ ตัวแปรเวลาในคอมพิวเตอร์จึงเดินจาก payoff กลับมาหาปัจจุบันของสัญญา

บนกริดราคาที่ห่างเท่ากัน เราประมาณ Delta ด้วย central difference และ Gamma ด้วย second difference

Delta บอกความชันของราคา Option ต่อราคาหุ้น ส่วน Gamma บอกว่าความชันเปลี่ยนเร็วเพียงใด ตัวอย่าง $\Delta S=5$ และค่า Option ที่ราคาสามจุดเป็น 7, 10, 14 จะได้ Delta≈0.7 และ Gamma≈0.04 ต่อดอลลาร์ ตัวเลขสามจุดนี้เป็นตัวอย่างสำหรับฝึกสูตร ไม่ใช่ผลลัพธ์กริดจริง

Forward difference ใช้ $(U_{i+1}-U_i)/\Delta S$ ส่วน backward difference ใช้ $(U_i-U_{i-1})/\Delta S$ เมื่อคำตอบเรียบพอ central difference ของ Delta และ Gamma มี truncation error อันดับ $O((\Delta S)^2)$ ขณะที่ forward/backward difference ของ Delta มีอันดับ $O(\Delta S)$ แต่ใกล้ขอบกริดหรือเมื่อ drift เด่นกว่า diffusion การเลือกด้านเดียวอาจเหมาะกว่าด้านเสถียรภาพ

สำหรับ Theta ต้องระวังเครื่องหมาย เพราะ $\partial_tV=-\partial_\tau U$

จาก Black–Scholes PDE สู่สูตรสามจุด

เขียน PDE ด้วยเวลา τ จะได้

แทนด้านซ้ายด้วย $(U_i^{k+1}-U_i^k)/\Delta\tau$ และใช้อันดับเวลา k สำหรับทุกพจน์ด้านขวา จะได้วิธี explicit: ค่าชั้นใหม่คำนวณจากค่าที่รู้แล้วสามจุด

เมื่อ $S_{\max}=400$, M=80, L=1,000 จะได้ ΔS=5 ดอลลาร์ และ Δτ=0.001 ปี ที่ S=100 หรือ i=20 สัมประสิทธิ์เป็น a=0.0077, b=0.98397, c=0.0083 ผลรวมเท่ากับ 0.99997 หรือ $1-r\Delta\tau$

ที่วันหมดอายุ Call มีค่า 0, 0, 5 ดอลลาร์บน S=95,100,105 ดังนั้นการเดินกลับจาก expiry หนึ่ง step ให้ค่าที่ S=100 เป็น $0.0077(0)+0.98397(0)+0.0083(5)=0.0415$ ดอลลาร์ แล้วใช้ทั้งแถวที่คำนวณได้ทำ step ถัดไป

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

สูตรนี้คล้ายการคิดย้อนกลับบนต้นไม้ แต่ a,b,c ไม่ใช่ probabilities ที่รวมได้หนึ่ง โดยตรง มีส่วนคิดลดรวมอยู่ด้วย และบางการตั้งค่าทำให้สัมประสิทธิ์ติดลบได้

ต้องรู้ payoff และค่าที่ขอบก่อนเริ่ม

ที่ τ=0 เรากำหนดทุกจุดเป็น payoff

นี่คือ terminal condition ในเวลา t และเป็น initial condition ในเวลา τ ส่วน boundary conditions คือค่าที่ปลายกริดราคาซึ่งต้องกำหนด ทุกชั้นเวลา

หุ้นที่เริ่มเป็นศูนย์จะคงเป็นศูนย์ใน GBM: Call จึงไม่มี payoff ส่วน Put จ่าย K เมื่อหมดอายุและมีค่าปัจจุบัน $Ke^{-r\tau}$ สำหรับ Call ที่ S สูงมาก มูลค่าเข้าใกล้หุ้นหนึ่งหน่วยลบมูลค่าปัจจุบันของ K

Smax เป็นการแทนช่วงราคาอนันต์ด้วยขอบจำกัด จึงเกิด domain truncation error ได้ หากขยับ Smax ออกแล้วราคาที่สนใจเปลี่ยนมาก แสดงว่าขอบเดิมใกล้เกินไป แต่ต้องคง ΔS ให้ใกล้เดิมด้วย มิฉะนั้นเราจะเปลี่ยนทั้งตำแหน่งขอบและความละเอียดพร้อมกัน

เอกสารยังเสนอขอบบนแบบ $\Gamma\approx0$ หรือ $U_M^k=2U_{M-1}^k-U_{M-2}^k$ สำหรับ payoff ที่เกือบเป็นเส้นตรงเมื่อ S สูง เป็นทางเลือกแบบ linear extrapolation ซึ่งต้องตรวจความเหมาะสมแยกต่างหาก ห้องทดลองนี้ใช้ค่าขอบ Call/Put ในตารางโดยตรง

กริดละเอียดขึ้นอาจทำให้วิธี explicit ใช้ไม่ได้

การใช้ time step ใหญ่เกินไปอาจขยายความผิดพลาดจนคำตอบสั่นหรือระเบิด สำหรับตัวอย่างที่ r≥0 เราใช้เกณฑ์เพียงพอที่ตรวจง่าย: สัมประสิทธิ์ a,b,c ของทุกจุดภายในต้องไม่ติดลบ ทำให้แต่ละ step ไม่ขยายค่าสูงสุดจากการถ่วงน้ำหนักภายในกริด เพราะผลรวมไม่เกินหนึ่ง

เงื่อนไขของ b ให้ข้อจำกัด

ที่ค่าตั้งต้น M=80 ต้องมี L อย่างน้อย 250 step เมื่อ T=1 ปี ดังนั้น 1,000 step ผ่านเงื่อนไขนี้ แต่ 100 step ไม่ผ่าน หากเพิ่ม M ประมาณสองเท่า time step ที่ยอมให้ใช้จะเล็กลงประมาณสี่เท่า งานคำนวณกริดหนึ่งมิติจึงเพิ่มเร็วกว่าจำนวนจุดราคาเพียงอย่างเดียว

ยังมีอีกเงื่อนไขที่การลด time step แก้ไม่ได้: $a_i\geq0$ ต้องการ $\sigma^2i\geq r$ ตัวอย่าง r=5%, σ=20%, i=1 จะได้ $\sigma^2i=0.04<0.05$ ทำให้ a ติดลบสำหรับทุก Δτ>0 ทางเลือกหนึ่งคือเปลี่ยนวิธีประมาณ drift หรือใช้กริด/วิธีที่เหมาะสม ไม่ใช่เพิ่ม L ไปเรื่อย ๆ

การไม่ผ่านเกณฑ์นี้ไม่ได้พิสูจน์ว่าทุกการรันจะระเบิด แต่ ไม่ผ่านเงื่อนไขรับรองที่ห้องทดลองเลือกใช้ จึงหยุดแสดงราคาและอธิบายเหตุผล เราไม่ต้องรอให้เห็นตัวเลขผิดปกติก่อนตรวจ scheme

เสถียรภาพและความแม่นยำก็เป็นคนละเรื่อง กริดที่ผ่านเกณฑ์อาจยังหยาบ ขอบราคาอาจใกล้เกินไป หรือ payoff อาจไม่เรียบที่ strike จึงต้องตรวจ convergence ต่อ

ทดลอง: เปลี่ยนกริดและตรวจราคาก่อนใช้

เริ่มจาก M=80, L=1,000 เทียบราคาและ Greeks กับ Black–Scholes จากนั้นลด L เป็น 100 เพื่อดูการตรวจเสถียรภาพ ลองเพิ่ม M โดยคง L แล้วใช้ปุ่มปรับจำนวน time step ให้ผ่านเกณฑ์ หากยังไม่ผ่าน ให้ดูเงื่อนไขด้าน drift ด้วย

ตรวจทั้งราคา Greeks และการลู่เข้า

เมื่อมีคำตอบบนกริด เราคำนวณ Delta และ Gamma จากจุดข้างเคียงได้ ส่วน Theta ที่วันนี้ประมาณจากสองชั้นเวลาสุดท้ายด้วย $(U_i^{L-1}-U_i^L)/\Delta\tau$ ซึ่งมีหน่วยดอลลาร์ต่อปี ห้องทดลองใช้ linear interpolation หาก S ไม่ตรงจุดกริด ส่วน Greeks ที่ขอบไม่แสดง เพราะสูตรกลางต้องมีเพื่อนบ้านสองด้าน

Gamma หักล้างค่าราคาที่อยู่ใกล้กันแล้วหารด้วย $(\Delta S)^2$ จึงไวต่อความผิดพลาดของกริด ราคาใกล้ benchmark ไม่ได้แปลว่า Greeks แม่นในระดับเดียวกัน โดยเฉพาะใกล้ strike และใกล้หมดอายุ

การตรวจที่มีประโยชน์คือ

1. ตรวจโจทย์ที่รู้คำตอบ: เทียบ Call/Put และ Greeks กับสูตร Black–Scholes
2. ลด Δτ โดยคง ΔS: ดูผลจากการแบ่งเวลา ภายใต้กริดที่ผ่านเกณฑ์
3. ลด ΔS พร้อมรักษาเสถียรภาพ: ดูว่าคำตอบเข้าใกล้กันหรือไม่
4. ขยาย Smax โดยคง ΔS: ตรวจผลของขอบจำกัด
5. ตรวจสมบัติของราคา: เช่น ไม่ติดลบ ขอบล่างตาม no-arbitrage และ $C-P=S-Ke^{-rT}$ โดยคำนึงถึง numerical tolerance

การประมาณ PDE มี formal truncation error $O(\Delta\tau+(\Delta S)^2)$ เมื่อคำตอบเรียบพอ แต่ payoff ของ Call/Put มีมุมที่ strike และเงื่อนไขขอบเป็นอีกแหล่ง error จึงไม่ควรสรุป global accuracy จากอันดับของสูตรเพียงอย่างเดียว

เลือกวิธีให้เข้ากับสัญญา

สำหรับ American Option แนวคิดบนกริดคือเปรียบเทียบค่าถือสัญญาต่อกับค่าที่ได้จากการใช้สิทธิทันที: $U_i^{k+1}=\max(U_{i,\mathrm{continue}}^{k+1},g(S_i))$ ต้องกำหนดขอบและ scheme ให้สอดคล้องกับปัญหา American ด้วย ห้องทดลองในหน้านี้และ Notebook ใช้ European Call/Put เท่านั้น

วิธี explicit เริ่มเขียนและตรวจได้ง่าย แต่ติดข้อจำกัด time step วิธี implicit และ Crank–Nicolson เป็นหัวข้อต่อยอดที่เปลี่ยนการแก้ชั้นเวลา และมีข้อพิจารณาด้านความแม่นยำ/การสั่นใกล้ payoff ของตนเอง

ลองตอบก่อนเปิด Notebook

1. Monte Carlo ให้ SE=0.12 ดอลลาร์จาก 10,000 รอบ ต้องใช้ประมาณกี่รอบให้ SE เหลือ 0.03 ดอลลาร์?

2. ถ้าเพิ่มจำนวนช่วงราคา M เป็นสองเท่า แต่คง L เดิม ทำไมกริด explicit ที่เคยผ่านเกณฑ์อาจไม่ผ่าน?

3. ทำไมการใช้ exact GBM ที่วันต้นและวันปลายสัญญายังไม่พอสำหรับ continuous barrier?

4. ที่ S=0 ทำไม European Put มีค่า $Ke^{-r\tau}$ แทน K และเหตุใดเงื่อนไขนี้จึงไม่ใช้กับ American Put แบบตรง ๆ?

ดาวน์โหลด Python Notebook เพื่อรัน Monte Carlo ตรวจ SE และเปรียบเทียบ finite difference หลายกริดกับ Black–Scholes โค้ดใช้ Python standard library และฝังภาพประกอบไว้ในไฟล์

แหล่งที่มาและขอบเขต

เรียบเรียงจากเอกสาร Introduction to Numerical Methods ของหลักสูตร Certificate in Quantitative Finance ในไฟล์ JA253.4 Notes.pdf ที่ผู้ใช้ให้มา 62 หน้า: risk-neutral Monte Carlo หน้า 4–17, กริดและอนุพันธ์หน้า 18–34, payoff และ explicit scheme หน้า 35–48, boundary conditions หน้า 49–58 และการเลือกวิธีหน้า 59–62

เนื้อหานี้อธิบายใหม่เป็นภาษาไทย เพิ่มตัวอย่างคำนวณ ห้องทดลองและการตรวจเสถียรภาพ ไม่ได้เผยแพร่ไฟล์ PDF ต้นฉบับหรือภาพหน้าสไลด์ จุดขยายความคือความต่างของ weak/strong error, exact update กับการสังเกตเส้นทาง, ช่วงความเชื่อมั่นของ Monte Carlo และข้อจำกัดสัมประสิทธิ์ central difference

เอกสารประกอบสำหรับตรวจสูตร: Mike Giles, University of Oxford, Monte Carlo Lecture 1 สำหรับ terminal GBM และ Box–Muller และ Lecture 9 สำหรับ strong/weak convergence

ราคาและกราฟทุกชุดเป็นตัวอย่างสมมติภายใต้แบบจำลอง ไม่ใช่ market quote หรือการสอบเทียบกับตลาด รายละเอียดวิธีคำนวณและแหล่งที่มาอยู่ใน provenance