Bài 15: Phương pháp sai phân hữu hạn (Finite Difference Method) trong định giá quyền chọn

Chúng ta đã học cây nhị phân và mô phỏng Monte Carlo, hai trong ba phương pháp số không thể thiếu trong Quant Finance hiện đại. Bài viết này giới thiệu phương pháp thứ ba, cũng lâu đời và quan trọng không kém, phương pháp sai phân hữu hạn (Finite Difference Method, FDM).

Bạn đọc có thể tham khảo thêm cuốn sách “Finite Difference Methods in Financial Engineering” Daniel J. Duffy (2006), tác giả trình bày rất chi tiết và gần như đề cập tới mọi khía cạnh của phương pháp này [1].

Trong bài viết này:

• 1. Từ Binomial Tree, Monte Carlo đến Finite Difference Method
• 2. Nền tảng lý thuyết của Finite Difference Method
  • 2.1. Black-Scholes PDE
  • 2.2. Định lý Feynman-Kac: Cầu nối giữa PDE và kỳ vọng
• 3. Finite Difference Method trong định giá quyền chọn
  • 3.1. Xây dựng lưới sai phân hữu hạn
  • 3.2. Điều kiện cuối và điều kiện biên
  • 3.3. Xấp xỉ đạo hàm: Ba phương pháp sai phân
• 4. Explicit Method (FTCS)
  • 4.1. Xây dựng công thức Explicit từ PDE
  • 4.2. Thực hành tính toán thủ công
  • 4.3. Thực hành Python
  • 4.4. Điều kiện ổn định von Neumann
• 5. Implicit Method (BTCS)
  • 5.1. Xây dựng công thức Implicit từ PDE
  • 5.2. Hệ phương trình ba đường chéo và thuật toán Thomas
  • 5.3. Thực hành tính toán thủ công
  • 5.4. Thực hành Python
• 6. Crank-Nicolson Method
  • 6.1. Ý tưởng: Trung bình cộng của hai thế giới
  • 6.2. Độ chính xác bậc hai
  • 6.3. Hiện tượng dao động
• 7. So sánh ba phương pháp
• 8. Tóm tắt và thảo luận
• Appendix

1. Từ Binomial Tree, Monte Carlo đến Finite Difference Method

Hãy hình dung cả ba phương pháp số đều đang đi tìm lời giải cho cùng một câu hỏi: “Quyền chọn này đáng giá bao nhiêu?”, nhưng theo ba góc nhìn hoàn toàn trái ngược:

• Cây nhị phân (Bài 6-9) chia thời gian thành từng bước rời rạc, tại mỗi bước tung một “đồng xu” nhằm mô phỏng giá tài sản tăng lên hoặc hạ xuống, rồi kiên nhẫn tính ngược từ tương lai về hiện tại.
• Monte Carlo (Bài 10-14) mô phỏng hàng chục nghìn kịch bản ngẫu nhiên của đường đi giá tài sản tiến tới tương lai, rồi gom tất cả lại để lấy một con số trung bình.
• FDM tư duy hoàn toàn khác biệt: giải trực tiếp phương trình vi phân riêng phần (Partial Differential Equation, PDE) mô tả giá quyền chọn, bằng cách rời rạc hóa cả trục giá tài sản lẫn trục thời gian thành một tấm lưới (grid), rồi tính giá trị tại từng nút của tấm lưới đó.

Black-Scholes (1973) đã cho công thức đóng để định giá European Option, vậy tại sao vẫn cần giải PDE? Câu trả lời giống hệt lý do chúng ta cần Monte Carlo: công thức đóng chỉ tồn tại cho một số ít trường hợp lý tưởng. Ngay khi điều kiện phức tạp hơn một chút, như khi định giá quyền chọn kiểu Mỹ, hay khi định giá với độ biến động thay đổi, PDE tương ứng không còn nghiệm giải tích nữa. Lúc này, FDM trở thành công cụ để giải gần đúng phương trình vi phân đó.

Trong ba phương pháp, FDM thường được coi là tiêu chuẩn công nghiệp để định giá các sản phẩm có quyền thực hiện sớm hoặc phụ thuộc vào một vài biến ngẫu nhiên (như lãi suất và giá cổ phiếu cùng lúc), nhờ tốc độ và độ chính xác vượt trội trong không gian ít chiều [1]. Nhưng cũng như mọi công cụ khác, FDM có những giới hạn riêng, chúng ta sẽ bàn kỹ ở Bài 17.

2. Nền tảng lý thuyết của Finite Difference Method

2.1. Black-Scholes PDE

Nhớ lại Bài 6. Binomial Model, ta xây dựng một danh mục mô phỏng gồm $\Delta$ cổ phiếu và $B$ tiền mặt, sao cho giá trị danh mục này luôn khớp với giá trị quyền chọn ở mọi kịch bản. Điều then chốt là danh mục này phải được tái cân bằng liên tục theo thời gian để đảm bảo luôn trung hòa rủi ro.

Bây giờ, hãy tưởng tượng ta thu nhỏ bước thời gian $\Delta t \to 0$ trong mô hình nhị phân đó, giống như cách chúng ta đã biến bước đi ngẫu nhiên thành chuyển động Brown ở Bài 4. Brownian Motion. Các nhà toán học dự đoán được phương trình sai phân rời rạc sẽ hội tụ về một phương trình vi phân liên tục. Đó chính xác là những gì Fischer Black, Myron Scholes, và Robert Merton đã làm vào năm 1973, dùng lập luận danh mục mô phỏng và phòng hộ liên tục theo thời gian để suy ra một phương trình vi phân riêng phần mà giá quyền chọn $V(S,t)$ theo giá tài sản cơ sở $S$ và thời gian $t$, bắt buộc phải thỏa mãn [1]:

Trong đó:

• $r$ lãi suất phi rủi ro.
• $\frac{\partial V}{\partial t}$ đạo hàm riêng bậc nhất của giá quyền chọn $V$ theo thời gian $t$.
• $\frac{\partial V}{\partial S}$, $\frac{\partial^2 V}{\partial S^2}$ đạo hàm riêng bậc nhất và bậc hai của giá quyền chọn $V$ theo giá tài sản $S$.

Đây chính là phương trình vi phân riêng phần Black-Scholes, gọi tắt là Black-Scholes PDE. Lưu ý cần phân biệt phương trình này với công thức Black-Scholes trong định giá quyền chọn. Đạo hàm trong phương trình là đạo hàm riêng phần $\partial$ vì $V$ phụ thuộc vào hai biến độc lập, giá cổ phiếu $S$ và thời gian $t$; chứ không phải đạo hàm toàn phần $d$ chỉ phụ thuộc $t$, ví dụ như phương trình vi phân ngẫu nhiên GBM mà ta đã học $dS_t = \mu S_t dt + \sigma S_t dW_t$.

Xem chứng minh tại Phụ lục (A1), hoặc bạn đọc có thể tạm thời bỏ qua vì phần chứng minh cần sử dụng đến Bổ đề Itô chưa học. Về mặt trực giác, phương trình (15.1) chính là phiên bản liên tục của hệ phương trình mà ta đã giải ở Bài 6 \(C_u = \Delta uS + r_0 B\) và \(C_d = \Delta dS + r_0 B\), khi số bước của cây tiến tới vô cùng.

Ý nghĩa của từng số hạng

Bằng cách chuyển $rV$ sang vế phải, phương trình (15.1) có thể viết lại như sau:

Để hiểu phương trình (15.2) một cách tự nhiên nhất, hãy tưởng tượng bạn bỏ ra một số tiền để mua quyền chọn và thực hiện một chiến lược phòng hộ hoàn hảo (vừa mua quyền chọn, vừa bán khống cổ phiếu để triệt tiêu mọi rủi ro biến động giá). Ngược lại, nếu đem số tiền này gửi ngân hàng, mỗi năm bạn chắc chắn thu về khoản tiền lãi với lãi suất phi rủi ro $r$. Theo nguyên lý no-arbitrage, để không ai bị thiệt hay thu lời bất chính, thu nhập từ danh mục quyền chọn đã phòng hộ phải bằng khoản tiền lãi gửi ngân hàng.

Vế phải là tiền lãi ngân hàng:

• $rV$: Tiền lãi từ số tiền $V$ nếu không mua quyền chọn.

Vế trái là thu nhập thuần từ toàn bộ danh mục quyền chọn, bao gồm:

• $\dfrac{\partial V}{\partial t}$: Mất mát từ hao mòn thời gian. Đây chính là Theta trong ngôn ngữ Greeks, đại diện cho tốc độ giá quyền chọn giảm dần theo thời gian, ngay cả khi giá cổ phiếu đứng yên.
• $\dfrac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}$: Lợi nhuận thu được từ sự biến động của giá cổ phiếu, được đo bằng Gamma ($\Gamma = \frac{\partial^2 V}{\partial S^2}$) khuếch đại bởi độ biến động $\sigma$. Nếu Delta ($\Delta = \frac{\partial V}{\partial S}$) là tốc độ biến đổi của giá quyền chọn theo giá cổ phiếu, thì Gamma chính là gia tốc, đo lường tốc độ thay đổi của chính Delta. Biến động thị trường ($\sigma$) càng lớn, phần lợi nhuận từ Gamma này tạo ra càng nhiều. Nếu thị trường hoàn toàn đứng yên ($\sigma = 0$), số hạng này biến mất, và PDE quay về một phương trình vi phân đơn giản hơn nhiều.
• $rS\dfrac{\partial V}{\partial S}$: Để thực hiện phòng hộ, ta cần vay/gửi số tiền $S\frac{\partial V}{\partial S}$ để mua/bán cổ phiếu cơ sở. Đây chính là tiền lãi phải trả/được nhận từ khoản tiền phòng hộ này.

Như vậy về mặt trực giác, phương trình Black-Scholes thực chất là một phương trình cân bằng tài chính: thu nhập từ hai chiến lược đầu tư (hai vế của phương trình) phải bằng nhau.

2.2. Định lý Feynman-Kac: Cầu nối giữa PDE và kỳ vọng

Nhớ lại Bài 10. Monte Carlo Simulation, chúng ta đã làm quen với cách định giá quyền chọn bằng Monte Carlo, nơi giá quyền chọn là kỳ vọng dưới độ đo trung hòa rủi ro \(V = e^{-rT}\,\mathbb{E}^{\mathbb{Q}}[f(S_T)]\). Bây giờ Black-Scholes lại cho rằng $V$ cũng là nghiệm giải tích của phương trình (15.1). Đây có phải hai bài toán khác nhau tình cờ cho cùng đáp số, hay có một mối liên hệ sâu xa hơn giữa xác suất và giải tích?

Cơ sở toán học cho cầu nối này là Định lý Feynman-Kac. Định lý này phát biểu như sau:

• Xét quá trình ngẫu nhiên $S_t$ tuân theo một SDE có dạng $dS_t = \mu(S_t,t)\,dt + \sigma(S_t,t)\,dW_t$
• Cho $f(x)$ là hàm payoff, $V(x, t)$ là kỳ vọng có điều kiện \(V(S,t) = e^{-r(T-t)}\,\mathbb{E}[f(S_T) \mid S_t = S]\)
• Với điều kiện biên tại thời điểm cuối $V(S,T) = f(S)$
• Khi đó, $V(x, t)$ là nghiệm của PDE sau:

Khi áp dụng vào tài chính, tài sản $S_t$ được mô hình hóa tuân theo GBM dưới độ đo trung hòa rủi ro $\mathbb{Q}$, nên có thể thay $\mu(S,t) = rS$ và $\sigma(S,t) = \sigma S$. Bạn sẽ nhận ra ngay phương trình (15.3) chính là Black-Scholes PDE (15.1) ở phần trước. Phần chứng minh cho định lý, áp dụng cho quá trình $e^{-rt}V(S_t,t)$, được trình bày ở Phụ lục (A2) [3].

Để hiểu định lý Feynman-Kac một cách trực quan nhất mà không bị ngợp bởi các ký hiệu toán học, hãy hình dung định lý này như một cầu nối giữa hai cách tiếp cận khác nhau nhằm xác định giá trị hợp lý của hợp đồng quyền chọn. Hai cách tiếp cận như sau:

• Cách 1: Mô phỏng hàng triệu con đường ngẫu nhiên từ $t$ đến ngày đáo hạn $T$ để thu được phân phối xác suất của $S_T$. Tính dòng tiền $f(S_T)$ với từng kịch bản tại $T$, lấy trung bình cộng, sau đó chiết khấu về hiện tại.
• Cách 2: Ta biết giá trị quyền chọn tại thời điểm cuối $T$ là $V(S, T) = f(S)$. Thiết lập PDE mô tả sự biến thiên của $V$ nếu lùi lại một khoảng thời gian cực kỳ nhỏ $dt$ dựa trên drift và volatility, lặp lại và giải lùi từ $T$ về $t$ để tìm giá trị $V(S, t)$ ngay tại thời điểm hiện tại.

Nói cách khác, Feynman-Kac chứng minh rằng “kỳ vọng dưới độ đo rủi ro trung hòa” (cách tiếp cận xác suất của Monte Carlo) và “nghiệm PDE” (cách tiếp cận giải tích của FDM) là hai mặt của cùng một đồng xu. Dù bạn giải xuôi hay giải ngược, kết quả thu được là như nhau [3].

Đây là điểm khá thú vị, FDM không phải một công cụ độc lập, mà là góc nhìn thứ hai của cùng một bài toán mà cây nhị phân và Monte Carlo đã giải quyết bằng xác suất.

3. Finite Difference Method: Phương pháp sai phân hữu hạn

3.1. Xây dựng lưới sai phân hữu hạn

Để giải Black-Scholes PDE bằng phương pháp số, chúng ta phải thay thế các miền giá $S$ (có thể nhận bất kỳ giá trị dương nào) và miền thời gian $t$ (liên tục từ $0$ đến $T$), thành một tấm lưới (grid) hai trục với hữu hạn các điểm (grid point). Quá trình này, được gọi là rời rạc hóa, biến phương trình liên tục thành một hệ phương trình đại số mà máy tính có thể giải được.

• Trục giá cổ phiếu $S$: ta giới hạn miền giá trong khoảng $[0, S_{max}]$, với $S_{max}$ được chọn đủ lớn để coi như lớn vô cùng đối với bài toán, hay khiến giá quyền chọn OTM nhỏ tới mức không đáng kể. Chia khoảng này thành $N$ đoạn bằng nhau, mỗi đoạn có độ dài $\Delta S = S_{max}/N$. Các nút giá được đánh số $S_i = i \,\Delta S$, với $i = 0, 1, \ldots, N$.
• Trục thời gian $t$: chia khoảng $[0, T]$ thành $M$ đoạn bằng nhau, mỗi đoạn có độ dài $\Delta t = T/M$. Các nút thời gian được đánh số $t_j = j \, \Delta t$, với $j = 0, 1, \ldots, M$.

Kết quả là một tấm lưới hình chữ nhật gồm $(N+1) \times (M+1)$ nút, với mỗi nút $(i, j)$ đại diện cho một cặp giá trị $(S_i, t_j)$. Ký hiệu $V_{i,j}$ là giá trị xấp xỉ của giá quyền chọn tại nút lưới $(S_i, t_j)$, đây cũng chính là ẩn số ta cần tìm [2].

Việc lựa chọn các tham số của tấm lưới đòi hỏi những sự đánh đổi mà bạn cần cân nhắc kỹ lưỡng. Lưới càng mịn ($M$ và $N$ càng lớn), ước lượng càng chính xác nhưng lại đòi hỏi nhiều khối lượng tính toán và bộ nhớ hơn. Lưới phải đủ mịn để phản ánh chính xác hành vi của quyền chọn ở vùng gần giá thực hiện – nơi giá trị biến đổi nhanh nhất – nhưng nếu mở rộng lưới quá xa vào các vùng deep OTM hoặc deep ITM sẽ gây lãng phí. Tỉ lệ giữa $\Delta t$ và $(\Delta S)^2$ đặc biệt quan trọng đối với tính ổn định của sơ đồ, chúng ta sẽ phân tích ở phần sau.

3.2. Điều kiện cuối và điều kiện biên

Để giải PDE trên một miền hữu hạn, ta cần biết giá trị của $V$ tại các cạnh của tấm lưới. Đây gọi là điều kiện cuối (terminal condition) và điều kiện biên (boundary condition), tương tự như cách một bài toán vật lý cần biết nhiệt độ ban đầu và nhiệt độ tại hai đầu thanh kim loại để giải phương trình truyền nhiệt.

Điều kiện cuối (tại $j = M$, tức $t = T$): Chính là payoff của quyền chọn tại ngày đáo hạn, giống hệt bước khởi tạo mà ta đã làm ở bài về cây nhị phân:

Điều kiện biên $S=0$ (tại $i=0$): Nếu giá cổ phiếu bằng 0, theo mô hình GBM, nó sẽ mãi mãi bằng 0. Với call, quyền chọn vô giá trị $V_{0,j} = 0$. Với put, người nắm giữ chắc chắn sẽ thực hiện quyền và nhận $K$ tại ngày đáo hạn, nên giá trị hiện tại của nó chỉ đơn giản là $K$ chiết khấu về thời điểm $t_j$:

Điều kiện biên $S = S_{max}$ (tại $i=N$): Khi giá cổ phiếu cực lớn so với giá thực hiện $K$, call gần như chắc chắn sẽ được thực hiện với giá trị tiệm cận $S - K$ chiết khấu; còn put gần như chắc chắn vô giá trị:

Với ba điều kiện này, chúng ta có thể lấp đầy giá trị ở ba cạnh của tấm lưới: cạnh phải (điều kiện cuối), cạnh trên và cạnh dưới (điều kiện biên). Bài toán quy về vấn đề cơ bản: làm sao để tính các nút bên trong tấm lưới (cột $j = M-1, M-2, \ldots, 0$) dựa trên những gì đã biết? Đây chính là lúc chúng ta cần công cụ toán học tiếp theo để xấp xỉ đạo hàm.

3.3. Xấp xỉ đạo hàm: Ba phương pháp sai phân

Trước khi đến với các công thức toán học, hãy tưởng tượng một tình huống rất đời thường. Bạn đang đứng trên một ngọn đồi, muốn biết độ dốc (độ nghiêng) chính xác tại vị trí bạn đang đứng, nhưng bạn không có công cụ đo góc, chỉ có một cây thước dây và một máy đo độ cao.

• Cách 1 (Forward): Bạn đo độ cao tại vị trí hiện tại, đi tới một bước về phía trước (ví dụ 1 mét), đo độ cao ở đó, rồi lấy hiệu số độ cao chia cho 1 mét.
• Cách 2 (Backward): Bạn đo độ cao tại vị trí hiện tại, lùi lại một bước phía sau, đo độ cao ở đó, rồi tính hiệu số và chia cho 1 mét.
• Cách 3 (Central): Bạn đi tới một bước phía trước, đo độ cao; rồi lùi lại một bước phía sau, đo độ cao; lấy hiệu số hai điểm đó chia cho 2 mét, mà không cần đo tại chính vị trí bạn đang đứng.

Bằng trực giác, ta có thể đoán: nếu ngọn đồi không phải là một mặt phẳng hoàn toàn mà có độ cong như trong thực tế, thì cách 3 sẽ chính xác hơn hai cách còn lại, vì nó cân bằng thông tin từ cả hai phía, thay vì thiên lệch về một bên duy nhất. Hãy thử với một ví dụ số cụ thể sau đây để thấy rõ mức độ sai số khác nhau. Cho hàm $f(x) = e^x$, ta biết đạo hàm chính xác tại $x=1$ là $f’(1) = e^1 \approx 2.718282$. Áp dụng ba công thức xấp xỉ trên với bước nhảy $\Delta x$ khác nhau:

Quan sát cột sai số: mỗi khi $\Delta x$ giảm 10 lần, sai số của Forward cũng chỉ giảm khoảng 10 lần, trong khi sai số của Central giảm tới gần 100 lần. Đây là hệ quả trực tiếp của bậc sai số cắt cụt trong khai triển Taylor, điều mà ta sẽ cùng tìm hiểu ngay sau đây.

3.3.1. Forward Difference

Nền tảng toán học của kỹ thuật sai phân hữu hạn chính là khai triển Taylor, công cụ mà ta đã dùng để xấp xỉ $\log(1+x) \approx x$ ở Ôn Tập i. Log Return. Khai triển Taylor của $f(x+h)$ quanh $x$ với bước nhảy $h$ là:

Chuyển vế $h f’(x)$ rồi chia hai vế cho $h$:

Khi $h \to 0$, số hạng \(-\frac{h^2}{6}f'''(x)\) là một đại lượng bậc $\mathcal{O}(h^2)$ sẽ lớn hơn rất nhiều số hạng bậc cao hơn $\mathcal{O}(h^4)$, do đó toàn bộ phần dư có thể được gom gọn lại thành $\mathcal{O}(h^2)$ ở dòng 2. Tuy nhiên, phần dư bậc 2 này cũng bị áp đảo hoàn toàn bởi số hạng \(-\frac{h}{2}f''(x)\), bậc $\mathcal{O}(h)$, nên sai số khi xấp xỉ từ dòng 2 sang dòng 3, còn gọi là sai số cắt cụt, có bậc $\mathcal{O}(h)$.

3.3.2. Backward Difference

Khai triển Taylor của $f(x-h)$ quanh $x$ với bước nhảy nhỏ $h$ là:

Tương tự, chuyển vế $h f’(x)$ rồi chia hai vế cho $h$:

Như vậy, Backward Difference cũng có sai số cắt cụt bậc $\mathcal{O}(h)$.

3.3.3. Central Difference

Lấy (15.7) trừ đi (15.9):

Để ý các số hạng bậc chẵn \(h^2 f''(x)\) và $\mathcal{O}(h^4)$ đã bị triệt tiêu hoàn toàn khi trừ hai khai triển cho nhau. Chuyển vế $2h f’(x)$ rồi rồi chia hai vế cho $2h$, ta có:

Số hạng bị bỏ qua bây giờ là \(\frac{h^2}{6}f'''(x)\), có bậc $\mathcal{O}(h^2)$, chính xác hơn hẳn so với Forward/Backward. Đây là lý do vì sao Central Difference luôn được ưu tiên sử dụng khi xấp xỉ đạo hàm bậc nhất.

Black-Scholes PDE cần cả $\partial^2 V/\partial S^2$. Cộng (15.7) với (15.9), các số hạng bậc lẻ $h f’(x)$ và \(h^3 f'''(x)\) triệt tiêu lẫn nhau, ta có:

Chuyển vế \(h^2 f''(x)\) rồi chia hai vế cho $h^2$, ta có:

Như vậy, Central Difference cho đạo hàm bậc hai cũng có sai số bậc $\mathcal{O}(h^2)$.

3.3.4. Truncation Error

Thuật ngữ sai số cắt cụt (truncation error) chỉ đơn giản là: khi ta cắt khai triển Taylor vô hạn thành một xấp xỉ hữu hạn, phần bị cắt bỏ chính là sai số. Bảng dưới đây tóm tắt lại bậc sai số của ba công thức, thước đo quan trọng nhất khi đánh giá chất lượng một phương pháp số:

Đây chính là bộ công cụ chúng ta cần để rời rạc hóa Black-Scholes PDE. Bây giờ ta chỉ cần thay các đạo hàm $\partial V/\partial t$, $\partial V/\partial S$, $\partial^2 V/\partial S^2$ trong phương trình (15.1) bằng các công thức sai phân tương ứng trên tấm lưới đã xây dựng ở Phần 3.1.

4. Explicit Method (FTCS)

4.1. Xây dựng công thức Explicit từ PDE

Phương pháp đầu tiên và đơn giản nhất là Explicit Method, dựa trên sơ đồ FTCS (Forward-Time Central-Space) như Hình 3. Ta dùng Forward Difference để tính đạo hàm theo thời gian, và Central Difference để tính đạo hàm theo giá cổ phiếu.

Nhớ lại quy ước ở phần 3.1, ta cần lùi từ $j=M$ (ngày đáo hạn) về $j=0$ (hiện tại), giá trị quyền chọn $V_{i,j}$ tại cột $j$ chưa biết (số hạng màu đỏ), được ước lượng dựa trên các nút $V_{i-1,j+1}$, $V_{i,j+1}$, $V_{i+1,j+1}$ tại cột $j+1$, cột đã biết giá trị từ bước trước. Tại nút $(i, j+1)$, ta xấp xỉ các đạo hàm như sau:

Thay các xấp xỉ sai phân hữu hạn ở trên vào Black-Scholes PDE (15.1), thực hiện phép tính tại nút $(i, j+1)$ với $S_i = i \Delta S$, rồi giải $V_{i,j}$ (ẩn số duy nhất) xem Phụ lục (A3), ta được:

Cấu trúc của phương trình này cho thấy giá trị quyền chọn tại một nút bất kỳ bên trong lưới là một tổ hợp tuyến tính có trọng số của của ba nút lân cận tại bước thời gian tiếp theo – những giá trị mà chúng ta đã biết từ điều kiện biên hoặc từ các bước tính toán trước đó. Mỗi nút được tính độc lập, hoàn toàn tường minh (explicit), đó chính là nguồn gốc cái tên của phương pháp này [2].

Về mặt trực giác, tổng có trọng số này tương tự như một phép tính giá trị kỳ vọng: $a_i$ và $c_i$ đại diện cho các trọng số (xác suất giả) của việc giá tài sản giảm hoặc tăng, trong khi $b_i$ là trọng số của việc giá tài sản giữ nguyên, tất cả đều đã được điều chỉnh theo chiết khấu lãi suất phi rủi ro. Cách nhìn dưới dạng xác suất này giúp ta hình dung được sự liên kết giữa phương pháp sai phân hữu hạn với mô hình cây nhị phân đã học. Về bản chất, sơ đồ FTCS thực thi cùng một logic nhưng trên một lưới mịn hơn.

4.2. Thực hành tính toán thủ công

Hãy định giá một quyền chọn European Put với bộ tham số quen thuộc: $S_0 = 100$, $K = 100$, $T = 1$, $r = 5\%$, $\sigma = 20\%$. Giá Black-Scholes tham chiếu là $P_{BS} \approx 5.5735$.

Để dễ tính tay, ta chọn một lưới cực nhỏ: $N = 4$, và $M = 2$. Lưới giá có hai trục $S = {0, 50, 100, 150, 200}$ nên bước giá $\Delta S = 50$) và $t = {0, 0.5, 1}$ nên bước thời gian $\Delta t = 0.5$. Đây là một tấm lưới quá thô để cho kết quả chính xác, nhưng đủ nhỏ để bạn có thể tự tay kiểm chứng từng con số.

Trước hết, tính các hệ số $a_i, b_i, c_i$ theo (15.14):

Điều kiện cuối ($j=2$, $t=T=1$): Công thức (15.4) cho put như sau $V_{i,2} = \max(100 - S_i,\ 0)$. Áp dụng cho $i=0,1,2,3,4$, ta có:

Bước lùi 1 ($j=1$, $t=0.5$): Áp dụng công thức (15.14) cho $i=1,2,3$, và điều kiện biên (15.5), (15.6) cho $i=0,4$, ta có:

Bước lùi 2 ($j=0$, $t=0$): Lặp lại quy trình, dùng kết quả ở cột $j=1$ vừa tính:

Kết quả tại nút $(i=2, j=0)$ ứng với $S = 100$, $t = 0$, chính là giá quyền chọn ta cần $V_{2,0} \approx \mathbf{1.3838}$.

Giá quyền chọn này (1.3838) khác xa giá Black-Scholes (5.5735), điều này hoàn toàn dễ hiểu vì lưới quá thô, giống hệt trải nghiệm ở bài trước với cây nhị phân 1-2 bước. Mục đích của ví dụ tính tay không phải để có kết quả chính xác, mà để bạn hiểu tường tận cách tính toán ngược từ điều kiện cuối về hiện tại.

4.3. Thực hành Python

Để thực thi Explicit FDM, ta chia nhỏ thời gian và giá tài sản thành một lưới các nút $(S, t)$ và tính ngược giá trị quyền chọn từ ngày đáo hạn về hiện tại. Thuật toán bao gồm các bước sau đây:

• Khởi tạo lưới và điều kiện cuối: Thiết lập các bước $\Delta S$, $\Delta t$, các hệ số $a_i, b_i, c_i$ và tính giá trị payoff tại ngày đáo hạn $T$ cho tất cả các nút giá $S_i$, ví dụ cho put f = np.maximum(K - S, 0.0).
• Tính các nút lưới bên trong: Tại mỗi bước lùi thời gian $j$, áp dụng công thức sai phân cho toàn bộ các nút bên trong ($i = 1 \dots N-1$) cùng lúc V_new[1:N] = a[1:N] * V[0:N-1] + b[1:N] * V[1:N] + c[1:N] * V[2:N+1].
• Cập nhật điều kiện biên: Gán trực tiếp giá trị lý thuyết tại hai biên $S_0 = 0$ và $S_N = S_{\max}$ theo từng loại quyền chọn, ví dụ cho put V_new[0], V_new[N] = K * np.exp(-r * tau), 0.0.
• Nội suy tuyến tính tại giá ban đầu: Lặp lại quá trình lùi thời gian cho đến $t=0$, sau đó tìm khoảng giá chứa $S_0$ (idx) và nội suy tuyến tính giữa hai nút lân cận để thu được giá quyền chọn chính xác: price = V[idx] + (S0 - S[idx]) / dS * (V[idx+1] - V[idx]).

Hãy cùng định giá quyền chọn bán với các tham số sau: $S_0=100$, $K=100$, $T=1$, $r=5\%$, $\sigma=20\%$, $S_{max}=300$, cố định $N=120$ (số bước giá), rồi thử với ba mức $M$ (số bước thời gian) khác nhau.

Bạn có thể chạy thử trên Google Colab.

import numpy as np

def explicit_fdm(S0, K, T, r, sigma, Smax, N, M, option_type='call'):
    dS, dt = Smax / N, T / M
    i = np.arange(0, N + 1)
    S = i * dS
    V = np.maximum(S - K, 0.0) if option_type == 'call' else np.maximum(K - S, 0.0)

    # Coefficients a_i, b_i, c_i
    a = 0.5 * dt * (sigma**2 * i**2 - r * i)
    b = 1 - dt * (sigma**2 * i**2 + r)
    c = 0.5 * dt * (sigma**2 * i**2 + r * i)

    for j in range(M, 0, -1):
        # Vectorized across all interior nodes at once
        V_new = V.copy()
        V_new[1:N] = a[1:N] * V[0:N-1] + b[1:N] * V[1:N] + c[1:N] * V[2:N+1]
        
        # Boundary condition
        tau = T - (j - 1) * dt
        V_new[0], V_new[N] = (0.0, Smax - K * np.exp(-r * tau)) if option_type == 'call' else (K * np.exp(-r * tau), 0.0)
        
        V = V_new

    # Linear interpolation to read the price at S0
    idx = int(S0 / dS)
    price = V[idx] + (S0 - S[idx]) / dS * (V[idx+1] - V[idx])
    return price

Chạy thử với bộ tham số ở trên với nhiều mức $M$ khác nhau, giữ nguyên $N=120$:

# Parameters
S0, K, T, r, sigma, Smax = 100.0, 100.0, 1.0, 0.05, 0.20, 300.0

# Pricing
for M in [20, 50, 100, 1000, 3000]:
    price = explicit_fdm(S0, K, T, r, sigma, Smax, 120, M, 'put')
    print(f"N={120} | M={M:5d} (dt={T/M:.5f}) | Explicit={price:>24,.2f}")

N=120 | M=   20 (dt=0.05000) | Explicit=  -56,806,074,026,675.11
N=120 | M=   50 (dt=0.02000) | Explicit= -319,074,587,115,928.62
N=120 | M=  100 (dt=0.01000) | Explicit=          -79,124,686.79
N=120 | M= 1000 (dt=0.00100) | Explicit=                    5.56
N=120 | M= 3000 (dt=0.00033) | Explicit=                    5.56

Chỉ với $\Delta t$ hơi lớn một chút ($M=20, 50, 100$), kết quả bùng nổ thành những con số âm khổng lồ, hoàn toàn phi lý cho một quyền chọn bán. Nhưng chỉ cần giảm $\Delta t$ xuống đủ nhỏ ($M=1000$), kết quả đột ngột trở nên hợp lý và hội tụ về giá Black-Scholes thật (5.5735). Đây chính là tính không ổn định của phương pháp Explicit FDM, sẽ được bàn ở phần sau.

4.4. Điều kiện ổn định von Neumann

Explicit FDM có một nhược điểm lớn: nếu bạn chọn $\Delta t$ quá lớn so với $\Delta S$, kết quả có thể bùng nổ cho ra những con số hoàn toàn vô nghĩa. Hãy hình dung một chiếc micro đặt quá gần loa: một tiếng động rất nhỏ ban đầu đi vào micro, được loa khuếch đại lên, quay lại micro, khuếch đại tiếp, và cứ thế lặp lại cho đến khi trở thành tiếng hú chói tai. Quay lại bài toán, khi lưới thời gian quá thô, sai số làm tròn ở một nút bị khuếch đại qua mỗi bước lùi thời gian, và nếu hệ số khuếch đại lớn hơn 1, sai số sẽ bùng nổ theo cấp số nhân [1].

Vậy làm sao chọn được $\Delta t$ và $\Delta S$ an toàn cho thuật toán? Câu trả lời đến từ phân tích tính ổn định của John von Neumann ra đời từ những năm 1940 khi ông nghiên cứu các phương trình sai phân cho bài toán khí động lực học. Phần chứng minh đầy đủ cần sử dụng phân tích Fourier tôi xin phép bỏ qua, kết quả cuối cùng cho Explicit Method áp dụng lên Black-Scholes PDE là [1]:

Với $N = \frac{S_{\max}}{\Delta S}$, bất đẳng thức trên có thể đơn giản hóa như sau:

Trực giác của điều kiện này nằm ngay trong hệ số $b_i = 1 - \Delta t(\sigma^2 i^2 + r)$ trong công thức (15.14). Để $b_i \geq 0$ tại nút khó nhất ($i=N$, ở biên trên của lưới, nơi $i^2$ lớn nhất), ta cần $\Delta t (\sigma^2 N^2 + r) \leq 1$. Nếu vi phạm điều kiện này, $b_i$ trở nên âm và đủ lớn về độ lớn tuyệt đối để khiến sai số bị khuếch đại qua mỗi vòng lặp thay vì suy giảm.

Hãy kiểm chứng lại với thực nghiệm ở trên với $N=120$, $\sigma=0.2$, $r=0.05$. Ngưỡng an toàn là:

Ngưỡng an toàn hoàn toàn khớp với kết quả thực nghiệm. Với $M=1000$ $\Delta t = 0.001 < 0.001736$, kết quả ổn định hội tụ về giá Black-Scholes; với $M=100$ $\Delta t = 0.01 > 0.001736$, vi phạm điều kiện nên kết quả bùng nổ thành con số hoàn toàn vô nghĩa.

Điều kiện (15.15) cho ta một kết luận quan trọng: $\Delta t$ phải nhỏ hơn một ngưỡng tỷ lệ nghịch với $\Delta S^2$. Nói cách khác, nếu muốn tăng độ chính xác lên gấp đôi bằng cách giảm bước giá $\Delta S$ đi một nửa $(\Delta S \to \frac{\Delta S}{2})$, thì $\Delta t$ phải giảm đi 4 lần chỉ để đảm bảo tính ổn định. Đây chính là cái giá phải trả của Explicit FDM, và là động lực chính cho các phương pháp Implicit ở phần sau.

5. Implicit Method (BTCS)

5.1. Xây dựng công thức Implicit từ PDE

Điều kiện ổn định ngặt nghèo của Explicit Method dẫn chúng ta đến phương pháp thứ hai Implicit Method dựa trên sơ đồ BTCS (Backward-Time, Central-Space).

Ý tưởng rất đơn giản nhưng mang lại hiệu quả to lớn: thay vì ước lượng đạo hàm theo $S$ tại cột đã biết ($j+1$), ta ước lượng chúng (các số hạng màu đỏ) tại chính cột đang cần tìm ($j$):

Khác biệt duy nhất với công thức Explicit (15.13) là các số hạng bên phải giờ mang chỉ số $j$ thay vì $j+1$. Nhưng chính sự thay đổi tưởng chừng nhỏ này lại thay đổi hoàn toàn bản chất bài toán. Thay (15.17) vào Black-Scholes PDE và sắp xếp lại, xem chi tiết ở Phụ lục (A4), ta được:

Sự khác biệt căn bản ở đây là: trong công thức Explicit (15.14), vế phải chỉ chứa một ẩn số $V_{i,j}$, nên ta có thể tính trực tiếp. Nhưng trong công thức Implicit (15.18), vế trái chứa ba ẩn số chưa biết cùng lúc $V_{i-1,j}$, $V_{i,j}$, $V_{i+1,j}$, vì chúng đều thuộc cột $j$ mà ta đang cần tìm, chỉ có giá trị tại cột $j+1$ là ta đã biết. Phương trình không cho ta cách tính trực tiếp từng nút, mà ngầm định (implicit) một mối quan hệ giữa nhiều ẩn số cùng lúc, đòi hỏi ta phải giải một hệ phương trình.

5.2. Ma trận ba đường chéo và thuật toán Thomas

Viết lại công thức (15.18) cho toàn bộ các nút bên trong lưới từ $i=1,\ldots,N-1$ tại một cột $j$ cố định, ta thu được một hệ $N-1$ phương trình tuyến tính với $N-1$ ẩn số. Biểu diện dưới dạng ma trận $\mathbf{A}\mathbf{\vec{x}} = \mathbf{\vec{d}}$:

Vector $\mathbf{\vec{x}}$ là vector cột chứa $N-1$ giá trị giá quyền chọn tại thời điểm $j$ cần giải. Vectơ $\mathbf{\vec{d}}$ chứa toàn bộ thông tin đã biết từ bước thời gian tương lai $j+1$ kết hợp với thông số từ hai biên $V_{0,j}$ và $V_{N,j}$.

Ma trận $\mathbf{A}$ có kích thước $N-1 \times N-1$, chứa toàn bộ các hệ số $\alpha_i, \beta_i, \gamma_i$ đứng trước các ẩn số $V_{i,j}$ tại thời điểm $j$. Ma trận này là ma trận ba đường chéo (tridiagonal matrix) với một cấu trúc đặc biệt: đường chéo chính $\beta_i$ tương ứng với nút trọng tâm $V_{i,j}$ đang xét, hai đường chéo liền kề $\alpha_i$, $\gamma_i$ tương ứng với nút biên dưới $V_{i-1,j}$ và biên trên $V_{i+1,j}$ lân cận, tất cả các phần tử khác đều bằng 0. Cấu trúc đặc biệt này cho phép ta giải hệ phương trình này cực kỳ hiệu quả, không cần dùng các phương pháp giải ma trận tổng quát thông thường rất tốn kém với độ phức tạp thuật toán $\mathcal{O}(N^3)$.

Thuật toán Thomas (TDMA – TriDiagonal Matrix Algorithm) khai thác đúng cấu trúc này để giải hệ phương trình trong thời gian $\mathcal{O}(N)$. Ý tưởng của thuật toán gồm hai giai đoạn:

• Giai đoạn tính xuôi đi từ nút đầu tiên đến nút cuối cùng, dùng phép khử Gauss để loại bỏ hệ số $\alpha_i$ trên đường chéo dưới khỏi mỗi phương trình, biến ma trận thành dạng tam giác trên đơn giản hơn.
• Giai đoạn tính ngược giải từ nút cuối cùng về nút đầu tiên, mỗi phương trình chỉ còn một ẩn số duy nhất.

Trong thực hành, các thư viện như scipy.linalg.solve_banded trong Python đã có sẵn phiên bản tối ưu của thuật toán, mà chúng ta sẽ dùng trực tiếp ở phần sau.

5.3. Thực hành tính toán thủ công

Để so sánh với Explicit Method, ta xây dựng tấm lưới với cùng tham số $S_{max} = 200$, $N = 4$ ($\Delta S = 50$), $M = 2$ ($\Delta t = 0.5$). Trước hết, tính các hệ số $\alpha_i, \beta_i, \gamma_i$ theo (15.18):

So sánh Bảng 5 và Bảng 4 ở phần trước, ba hệ số Implicit chỉ là phiên bản đối xứng của ba hệ số Explicit: $\alpha_i = -a_i$, $\gamma_i = -c_i$ và $\beta_i = 2 - b_i$, nhờ vậy ta không cần tính lại từ đầu. Với Explicit, mỗi nút của cột $j$ được tính trực tiếp từ cột $j+1$. Với Implicit, ba ẩn số $V_{1,j}, V_{2,j}, V_{3,j}$ trên cùng một cột ràng buộc lẫn nhau, buộc ta phải giải đồng thời cả hệ ba phương trình (15.19).

Điều kiện cuối ($j=2$, $t=T=1$) giống hệt phần Explicit:

Bước lùi 1 ($j=1$, $t=0.5$). Áp dụng công thức (15.18) cho $i=1,2,3$. Điều kiện biên cho $V_{4,1} = 0$ và $V_{0,1} = 100\,e^{-0.05 \times 0.5} \approx 97.5310$. Hai giá trị này đã biết nên có thể chuyển số hạng chứa chúng $\alpha_1\, V_{0,1}$, $\gamma_3\, V_{4,1}$ sang vế phải. Ta có hệ phương trình ba ẩn (màu đỏ) như sau:

Thay giá trị đã biết rồi giải hệ phương trình, ta có:

Một cách kiểm tra nhanh nghiệm của hệ phương trình là thay ngược vào phương trình đầu tiên. Kết quả khớp: $1.045 \times 47.6275 - 0.0225 \times 0.6482 = 49.7562 \approx 50 - 0.0025 \times 97.5310$.

Trong ví dụ này, ta có thể viết dưới dạng hệ phương trình và giải tay một cách nhanh chóng. Tuy nhiên việc giải tay sẽ bất khả thi nếu số lượng biến tăng lên đáng kể, lúc này bài toán cần được viết lại dưới dạng ma trận để tính toán thuận tiện hơn.

Bước lùi 2 ($j=0$, $t=0$). Ma trận hệ số không đổi, các hệ số $\alpha_i, \beta_i, \gamma_i$ không phụ thuộc thời gian. Điều kiện biên mới cho $V_{4,0} = 0$ và $V_{0,0} = 100\,e^{-0.05} \approx 95.1229$. Ta có hệ phương trình ba ẩn mới như sau:

Kết quả tại nút $(i=2, j=0)$ ứng với $S = 100$, $t = 0$, giá quyền chọn theo Implicit là $V_{2,0} \approx \mathbf{1.2070}$.

Cũng như Explicit (1.3838), Implicit cho kết quả (1.2070) sai lệch lớn với giá Black-Scholes (5.5735), vì lưới quá thô. Mục đích của ví dụ này là để bạn thấy mỗi bước lùi của Implicit phải giải một hệ phương trình thay vì tính từng nút.

5.4. Thực hành Python

Implicit Method bắt đầu bằng việc xây dựng lưới $(S, t)$, sau đó giải hệ phương trình tuyến tính ba đường chéo để ước tính giá trị quyền chọn từ tương lai lùi về hiện tại.

• Khởi tạo lưới và ma trận hệ số ba đường chéo: Thiết lập các khoảng chia $\Delta S$, $\Delta t$, tính payoff tại ngày đáo hạn $T$, đồng thời xây dựng ma trận A chứa các hệ số $\alpha_i, \beta_i, \gamma_i$.
• Thiết lập vector $\mathbf{\vec{d}}$ tích hợp điều kiện biên: Tại mỗi bước thời gian $j$, khởi tạo vector d từ các giá trị đã biết ở thời điểm $j+1$, tính giá trị lý thuyết ở hai biên $V_0, V_N$ và thêm vào hai đầu vector d[0] -= alpha[1] * V0 và d[-1] -= gamma[N-1] * VN.
• Giải hệ phương trình bằng Thuật toán Thomas: Áp dụng hàm solve_banded((1, 1), A, d) để giải hệ $\mathbf{A}\mathbf{\vec{x}} = \mathbf{\vec{d}}$, thu được nghiệm x cho các nút bên trong rồi ghép với hai biên để hoàn thiện cột giá trị: V = np.concatenate(([V0], x, [VN])).
• Nội suy tuyến tính tại giá ban đầu: Lặp lại quá trình lùi thời gian cho đến $t=0$, sau đó tìm khoảng giá chứa $S_0$ (idx) và nội suy tuyến tính giữa hai nút lân cận để thu được giá quyền chọn chính xác: price = V[idx] + (S0 - S[idx]) / dS * (V[idx+1] - V[idx]).

from scipy.linalg import solve_banded

def implicit_fdm(S0, K, T, r, sigma, Smax, N, M, option_type='call'):
    dS, dt = Smax / N, T / M
    i = np.arange(0, N + 1)
    S = i * dS
    V = np.maximum(S - K, 0.0) if option_type == 'call' else np.maximum(K - S, 0.0)

    # Coefficients alpha_i, beta_i, gamma_i
    alpha = 0.5 * dt * (r * i - sigma**2 * i**2)
    beta  = 1 + dt * (sigma**2 * i**2 + r)
    gamma = -0.5 * dt * (r * i + sigma**2 * i**2)

    # Tridiagonal matrix
    A = np.zeros((3, N - 1))
    A[0, 1:]  = gamma[1:N-1]   # superdiagonal
    A[1, :]   = beta[1:N]      # main diagonal
    A[2, :-1] = alpha[2:N]     # subdiagonal

    for j in range(M, 0, -1):
        # Boundary condition
        tau = T - (j - 1) * dt
        V0, VN = (0.0, Smax - K * np.exp(-r * tau)) if option_type == 'call' else (K * np.exp(-r * tau), 0.0)
        
        d = V[1:N].copy()
        d[0]  -= alpha[1] * V0
        d[-1] -= gamma[N-1] * VN
        x = solve_banded((1, 1), A, d) #Thomas algorithm
        V = np.concatenate(([V0], x, [VN]))

    idx = int(S0 / dS)
    price = V[idx] + (S0 - S[idx]) / dS * (V[idx+1] - V[idx])
    return price

Thử chạy lại với bộ tham số khiến thuật toán Explicit bùng nổ ở phần trên:

# Pricing
for M in [20, 50, 100]:
    ex_price = explicit_fdm(S0, K, T, r, sigma, Smax, 120, M, 'put')
    im_price = implicit_fdm(S0, K, T, r, sigma, Smax, 120, M, 'put')
    print(f"N={120} | M={M:3d} (dt={T/M:.5f}) | Explicit={ex_price:>24,.2f} | Implicit={im_price:>5,.2f}")

N=120 | M= 20 (dt=0.05000) | Explicit=  -56,806,074,026,675.11 | Implicit= 5.51
N=120 | M= 50 (dt=0.02000) | Explicit= -319,074,587,115,928.62 | Implicit= 5.54
N=120 | M=100 (dt=0.01000) | Explicit=          -79,124,686.79 | Implicit= 5.55

Ngay cả với $\Delta t$ rất lớn, Implicit Method vẫn cho ra những giá trị hợp lý dao động quanh giá Black-Scholes (5.5735), mà không hề gặp phải hiện tượng bùng nổ nào.

5.4. Ổn định vô điều kiện

Implicit Method phải vất vả giải một hệ phương trình tại mỗi bước thời gian để được gì? Câu trả lời là sự ổn định vô điều kiện. Hay nói cách khác, về mặt lý thuyết, ta có thể chọn $\Delta t$ lớn tùy ý mà không sợ kết quả bùng nổ như Explicit Method [1].

Để hiểu vì sao Implicit Method luôn đảm bảo tính ổn định, hãy tưởng tượng bạn đang dự báo thời tiết. Explicit giống như phép dự báo ngoại suy thô: nhìn vào xu hướng nhiệt độ hôm nay, rồi ngoại suy thẳng ra dự báo cho ngày mai mà không tham chiếu lại với các dự báo đó. Nếu xu hướng hôm nay biến động mạnh, phép ngoại suy này dễ dàng mắc lỗi, và sai lệch càng lúc càng xa thực tế qua từng ngày. Implicit giống như việc giải một hệ phương trình cân bằng: nhiệt độ ngày mai phải nhất quán với chính nó và với các điểm lân cận. Vì lời giải luôn được tìm sao cho tự nhất quán ở mọi điểm cùng lúc, nó không có cơ hội để sai lệch và tự khuếch đại qua từng bước.

Điều này không có nghĩa là Implicit Method hoàn hảo. Nhìn vào kết quả ở phần trước, khi $\Delta t$ lớn ($M=20$), kết quả Implicit (5.5114) vẫn có sai số đáng kể so với giá trị thật (5.5735), vì sai số cắt cụt của phương pháp vẫn còn đó. Ổn định chỉ đơn giản là thuật toán không đưa ra những con số vô nghĩa khi $\Delta t$ lớn, chứ không có nghĩa là kết quả sẽ chính xác. Cái giá phải trả cho sự ổn định này là chi phí tính toán: mỗi bước thời gian giờ đây cần giải một hệ phương trình dù đã được tối ưu bằng thuật toán Thomas, thay vì chỉ vài phép nhân đơn giản như Explicit Method.

6. Crank-Nicolson Method

6.1. Ý tưởng: Trung bình cộng của hai thế giới

Chúng ta vừa tìm hiểu hai cách tiếp cận: Explicit Method nhanh, đơn giản, nhưng có thể bùng nổ và Implicit Method ổn định vô điều kiện, nhưng cần giải hệ phương trình mỗi bước. Cả hai đều có sai số cắt cụt với bậc $\mathcal{O}(\Delta t)$. Liệu có cách nào kết hợp ưu điểm của cả hai, đồng thời cải thiện độ chính xác?

Phương pháp Crank-Nicolson được giới thiệu lần đầu bởi John Crank và Phyllis Nicolson năm 1947, một trong những kỹ thuật kinh điển và được sử dụng rộng rãi trong lĩnh vực giải tích số PDE. Phương pháp Crank-Nicolson thực hiện xấp xỉ đạo hàm theo không gian bằng cách lấy trung bình cộng giữa vế phải của Explicit Method và vế phải của Implicit Method [1] [4]:

Thay các xấp xỉ sai phân hữu hạn ở trên vào Black-Scholes PDE (15.1), thực hiện phép tính và nhóm các ẩn số chưa biết ở thời điểm $j$ sang vế trái, các giá trị đã biết ở thời điểm $j+1$ sang vế phải, xem chi tiết ở Phụ lục (A5), ta có:

trong đó $a_i, b_i, c_i$ chính là các hệ số đã định nghĩa ở phương trình (15.14).

Cấu trúc của phương trình này cho thấy tại mỗi bước lùi thời gian, vế trái tạo thành một hệ ma trận ba đường chéo chứa ba ẩn số $V_{i-1,j}, V_{i,j}, V_{i+1,j}$ cần giải đồng thời (tương tự sơ đồ Implicit), trong khi vế phải là một tổ hợp các giá trị đã biết tại cột $j+1$ (tương tự sơ đồ Explicit) và đóng vai trò là vector vế phải $\mathbf{d}$.

6.2. Độ chính xác bậc hai

Nhớ lại từ Phần 3.3, đạo hàm theo giá tài sản $\frac{\partial V}{\partial S}$ sử dụng Central Difference nên cho sai số cắt cụt bậc hai $\mathcal{O}(\Delta S^2)$, dù dùng ở Explicit hay Implicit. Nhưng riêng đạo hàm theo thời gian $\frac{\partial V}{\partial t}$, cả hai phương pháp dùng Forward/Backward Difference chỉ cho sai số bậc một $\mathcal{O}(\Delta t)$.

Phép lấy trung bình cộng trong công thức (15.20) thực chất tương đương với việc tính đạo hàm thời gian tại điểm giữa $t_{j+1/2} = (t_j + t_{j+1})/2$, như Central Difference. Crank-Nicolson vì vậy có thể đạt được sai số bậc hai $\mathcal{O}(\Delta t^2)$ theo thời gian thay vì chỉ $\mathcal{O}(\Delta t)$ như hai phương pháp trước, trong khi vẫn giữ nguyên sai số bậc hai $\mathcal{O}(\Delta S^2)$ theo giá [1].

6.3. Hiện tượng dao động

Crank-Nicolson có một nhược điểm cần lưu ý: khi điều kiện cuối là hàm không trơn, ví dụ như hàm payoff của call $\max(S-K,0)$, phương pháp này đôi khi tạo ra những dao động giả (spurious oscillations) quanh điểm gãy khúc trong vài bước thời gian đầu tiên, trước khi các dao động này tắt dần [1].

Hiện tượng này có nhiều tương đồng với hiệu ứng răng cưa mà chúng ta đã gặp ở Bài 9. Binomial Model Limitation, đều bắt nguồn từ việc một hàm không trơn bị ép vào một lưới rời rạc hữu hạn. Ở cây nhị phân, giải pháp là Black-Scholes Smoothing, thay payoff rời rạc tại bước ngay trước ngày đáo hạn bằng giá trị Black-Scholes mượt mà. Với Crank-Nicolson, một giải pháp có cùng ý tưởng là Rannacher smoothing (Rannacher, 1984), dùng vài bước Implicit thuần túy giúp làm mượt tại những bước đầu tiên gần đáo hạn, trước khi chuyển hẳn sang Crank-Nicolson cho các bước còn lại.

Không có phương pháp nào luôn luôn tốt nhất, Crank-Nicolson chính xác hơn với bậc sai số lý thuyết cao, nhưng đòi hỏi sự cẩn trọng nhất định khi payoff không trơn.

6.4. Thực hành Python

Khác biệt giữa phương pháp Crank-Nicolson và Implicit chính là cách thiết lập vector $\mathbf{\vec{d}}$. Tại mỗi bước thời gian $j$, vector d được tính bằng cách kết hợp trung bình trọng số của cả giá trị ở hiện tại lẫn tương lai ($j$ và $j+1$) d = a[1:N] * V[0:N-1] + (1 + b[1:N]) * V[1:N] + c[1:N] * V[2:N+1], sau đó cộng thêm các điều kiện biên tương ứng.

def crank_nicolson_fdm(S0, K, T, r, sigma, Smax, N, M, option_type='call'):
    dS, dt = Smax / N, T / M
    i = np.arange(0, N + 1)
    S = i * dS
    V = np.maximum(S - K, 0.0) if option_type == 'call' else np.maximum(K - S, 0.0)

    a = 0.25 * dt * (sigma**2 * i**2 - r * i) # a = a_explicit/2
    b = -0.5 * dt * (sigma**2 * i**2 + r)     # b = (b_explicit - 1)/2
    c = 0.25 * dt * (sigma**2 * i**2 + r * i) # c = c_explicit/2

    A = np.zeros((3, N - 1))
    A[0, 1:]  = -c[1:N-1]
    A[1, :]   = 1 - b[1:N]
    A[2, :-1] = -a[2:N]

    for j in range(M, 0, -1):
        # Boundary condition
        tau = T - (j - 1) * dt
        V0, VN = (0.0, Smax - K * np.exp(-r * tau)) if option_type == 'call' else (K * np.exp(-r * tau), 0.0)
        
        d = a[1:N] * V[0:N-1] + (1 + b[1:N]) * V[1:N] + c[1:N] * V[2:N+1]
        d[0]  += c[1] * V0
        d[-1] += a[N-1] * VN
        x = solve_banded((1, 1), A, d)
        V = np.concatenate(([V0], x, [VN]))

    idx = int(S0 / dS)
    price = V[idx] + (S0 - S[idx]) / dS * (V[idx+1] - V[idx])
    return price

Bây giờ, hãy so sánh cả ba phương pháp trên cùng một lưới sai phân:

# Pricing
bs_price = black_scholes(S0, K, T, r, sigma, 'put')

for N, M in [(60, 500), (60, 1000), (120, 1000), (120, 1500), (200, 1500)]:
    ex  = explicit_fdm(S0, K, T, r, sigma, Smax, N, M, 'put')
    im = implicit_fdm(S0, K, T, r, sigma, Smax, N, M, 'put')
    cn = crank_nicolson_fdm(S0, K, T, r, sigma, Smax, N, M, 'put')
    print(f"N={N:3d} M={M:4d} | Expl={ex:.4f} (err={abs(ex-bs_price):.4f}) | "
          f"Impl={im:.4f} (err={abs(im-bs_price):.4f}) | CN={cn:.4f} (err={abs(cn-bs_price):.4f})")

N= 60 M= 500 | Expl=5.5130 (err=0.0606) | Impl=5.5091 (err=0.0644) | CN=5.5110 (err=0.0625)
N= 60 M=1000 | Expl=5.5120 (err=0.0615) | Impl=5.5101 (err=0.0634) | CN=5.5110 (err=0.0625)
N=120 M=1000 | Expl=5.5590 (err=0.0145) | Impl=5.5571 (err=0.0164) | CN=5.5580 (err=0.0155)
N=120 M=1500 | Expl=5.5587 (err=0.0149) | Impl=5.5574 (err=0.0161) | CN=5.5580 (err=0.0155)
N=200 M=1500 | Expl=5.5780 (err=0.0045) | Impl=5.5768 (err=0.0032) | CN=5.5774 (err=0.0038)

Với lưới đủ mịn, cả ba phương pháp đều hội tụ rất gần giá trị thật, sai số chủ yếu đến từ việc chọn $S_{max}$ hữu hạn và độ mịn của $\Delta S$ ($N$ lớn), chứ không còn từ bậc hội tụ theo thời gian. Sự khác biệt chỉ thực sự bộc lộ khi ta ép $\Delta t$ lớn ($M$ nhỏ), như ta đã thấy ở Phần 4.3 và Phần 5.4.

for M in [20, 50, 100]:
    ex  = explicit_fdm(S0, K, T, r, sigma, Smax, 120, M, 'put')
    im = implicit_fdm(S0, K, T, r, sigma, Smax, 120, M, 'put')
    cn = crank_nicolson_fdm(S0, K, T, r, sigma, Smax, 120, M, 'put')
    print(f"M={M:3d} (dt={T/M:.5f}) | Explicit={ex:>24,.2f} | Implicit={im:.4f} | CN={cn:.4f}")

M= 20 (dt=0.05000) | Explicit=  -56,806,074,026,675.11 | Implicit=5.5114 | CN=5.5587
M= 50 (dt=0.02000) | Explicit= -319,074,587,115,928.62 | Implicit=5.5393 | CN=5.5581
M=100 (dt=0.01000) | Explicit=          -79,124,686.79 | Implicit=5.5487 | CN=5.5581

Đây chính là bằng chứng thực nghiệm rõ ràng nhất cho toàn bộ lý thuyết của bài viết này. Explicit sụp đổ hoàn toàn, trong khi Implicit và Crank-Nicolson vẫn cho ra những con số hợp lý, đặc biệt Crank-Nicolson nhờ độ chính xác bậc hai theo thời gian bám sát giá trị thật (5.5735) hơn hẳn Implicit ngay cả khi $\Delta t$ rất lớn.

7. So sánh ba phương pháp

Hãy tổng hợp lại những gì ba phương pháp mang lại, cả về mặt trực giác lẫn hiệu năng.

Trong thực tế công nghiệp, Crank-Nicolson (kết hợp với Rannacher smoothing cho vài bước đầu) thường là lựa chọn mặc định cho hầu hết bài toán định giá quyền chọn kiểu Âu và kiểu Mỹ đơn giản, nhờ sự cân bằng tốt nhất giữa tốc độ hội tụ và chi phí tính toán. Explicit Method hiếm khi được dùng trong sản xuất thực tế do ràng buộc ổn định quá ngặt nghèo, nhưng lại là công cụ giảng dạy tuyệt vời vì tính trực quan và mối liên hệ với cây tam phân.

Bạn có thể tự tay điều chỉnh Widget dưới đây để trực quan hóa toàn bộ tấm lưới sai phân, xem cách giá trị quyền chọn được tính ngược từ điều kiện cuối bên phải về hiện tại bên trái, cũng như quan sát trực tiếp hiện tượng bùng nổ của Explicit Method khi $\Delta t$ vượt ngưỡng ổn định.

Mỗi ô trên lưới bên trái là một nút $V(S,t)$ của bảng giá: trục tung là giá $S$, trục hoành là thời gian $t$ với ngày đáo hạn $T$ nằm ở bên phải và hiện tại $t=0$ nằm ở bên trái. Màu ô càng đậm thì giá trị quyền chọn càng cao, ô xám là cột chưa được tính. Đường thẳng đậm từ trên xuống là vị trí cột đang tính của phép tính lùi, từ phải qua trái. Đồ thị bên phải là lát cắt $V(S)$ tại đúng cột đang tính.

• Kéo $N, M$ xuống mức thấp nhất, bạn sẽ thấy lưới trở lên rất thô.
• Bấm Play hoặc kéo thanh Backward sweep: Cột con trỏ dịch chuyển, và lưới sai phân được tô dần từ phải sang trái thể hiện cho việc tính lùi từ ngày đáo hạn về hiện tại.
• Chọn Explicit, đặt $N = 60, M = 10$ rồi bấm Play: Metric báo $\Delta t = 0.00694 > \Delta t^* = 0.0025$ và hiện chữ đỏ UNSTABLE. Bạn sẽ thấy lưới không còn chuyển màu mượt mà, mà xuất hiện các nút lưới có màu đậm bất thường đại diện cho giá cổ phiếu lớn vô lý. Các ô màu này thường xuất hiện ở các hàng trên cùng (nơi $S$ lớn, $i$ lớn) rồi lan dần xuống dưới, vì hệ số $b_i = 1-\Delta t(\sigma^2 i^2+r)$ chứa $i^2$ sẽ bị khuếch đại mạnh. Đồ thị phía bên phải xuất hiện răng cưa, hoặc thậm chí có những đoạn trống trơn nơi giá bùng nổ.
• Giữ nguyên $N, M$, đổi sang Implicit hoặc Crank-Nicolson: Lưới trở lại bình thường, màu chuyển mượt mà. Đồ thị bên phải cũng khớp hoàn toàn với giá Black-Scholes. Như vậy nguyên nhân nằm ở sơ đồ chứ không phải ở bài toán.

8. Tóm tắt và thảo luận

Trong bài viết này, chúng ta đã đặt nền móng cho viên gạch thứ ba và cuối cùng của bộ ba phương pháp số trong Quant Finance.

• Ba phương pháp sai phân khác nhau đánh đổi giữa tốc độ, độ chính xác, và sự ổn định:
  • Explicit: nhanh, đơn giản, nhưng đòi hỏi điều kiện ổn định ngặt nghèo.
  • Implicit: ổn định vô điều kiện, đổi lại phải giải một hệ phương trình tại mỗi bước.
  • Crank-Nicolson: kết hợp cả hai, đạt độ chính xác bậc hai theo thời gian, là lựa chọn mặc định phổ biến nhất trong thực tế, dù cần cẩn trọng với hiện tượng dao động quanh payoff không trơn.

• Về mặt lý thuyết, định lý Feynman-Kac cho thấy cách tiếp cận theo Black-Scholes PDE và cách tiếp cận Monte Carlo bằng kỳ vọng dưới độ đo trung hòa rủi ro thực chất là hai mặt của cùng một bài toán. Khai triển Taylor là nền tảng của mọi công thức sai phân hữu hạn. Forward/Backward Difference cho sai số $\mathcal{O}(\Delta x)$, trong khi Central Difference cho sai số $\mathcal{O}(\Delta x^2)$ chính xác hơn hẳn.

• Thực hành Python

Ở bài này, chúng ta mới chỉ định giá European Option, nơi FDM không có ưu thế nổi trội nào so với công thức Black-Scholes. Bài tiếp theo chúng ta sẽ dùng công cụ này để chinh phục những sản phẩm phái sinh phức tạp hơn.

Appendix

A1. Chứng minh Black-Scholes PDE từ danh mục phòng hộ

Ta xây dựng một danh mục $\mathcal{\Pi}$ gồm bán khống 1 quyền chọn giá trị $V$ và mua $\Delta$ cổ phiếu:

Theo giả định, giá cổ phiếu tuân theo GBM dưới độ đo thực $dS = \mu Sdt + \sigma SdW_t$. Áp dụng Bổ đề Itô cho hàm $V(S_t, t)$, ta có:

Thay đổi giá trị của danh mục $\mathcal{\Pi}$ trong một khoảng thời gian $dt$ là:

Thay (15.A1.2) và SDE $dS = \mu Sdt + \sigma SdW_t$ vào (15.A1.3), ta có:

Nhóm các số hạng $dt$ và $dW_t$ riêng:

Đây chính là bước chọn $\Delta$ đã được giới thiệu trong phần danh mục mô phỏng ở Bài 6. Chọn $\Delta = \frac{\partial V}{\partial S}$ để triệt tiêu hoàn toàn số hạng ngẫu nhiên $dW_t$, đồng thời cũng triệt tiêu luôn số hạng $\mu S$ còn sót lại trong ngoặc $dt$, vì cả hai đều có chung thừa số $\left(\Delta - \frac{\partial V}{\partial S}\right)$. Phương trình (15.A1.3) lúc này trở thành:

Vì $d\Pi$ không còn chứa $dW_t$, danh mục $\mathcal{\Pi}$ trở thành phi rủi ro trong khoảng thời gian $dt$. Theo giả định no-arbitrage, một danh mục trung hòa rủi ro chỉ có thể sinh lời đúng bằng lãi suất phi rủi ro:

Cho (15.A1.4) bằng (15.A1.5):

Chuyển vế và sắp xếp lại:

Điều cần chứng minh, công thức (15.1).

Một điểm thú vị là lợi suất kỳ vọng của cổ phiếu $\mu$ đã hoàn toàn biến mất khỏi phương trình cuối cùng, đúng như ta đã học giá quyền chọn không phụ thuộc vào khẩu vị rủi ro của nhà đầu tư.

A2. Định lý Feynman-Kac: Chứng minh mối liên hệ giữa PDE và kỳ vọng

Ta chỉ chứng minh chiều quan trọng nhất của định lý cho mục đích của bài viết: nếu $V(S,t)$ thỏa mãn PDE (15.1) cùng điều kiện cuối $V(S,T)=f(S)$, thì $V(S,t) = e^{-r(T-t)}\,\mathbb{E}[f(S_T)\mid S_t=S]$.

Xét quá trình ngẫu nhiên $Y_t = e^{-rt}V(S_t, t)$. Áp dụng quy tắc tích của Bổ đề Itô cho tích của hai quá trình ngẫu nhiên $e^{-rt}$ và $V(S_t,t)$:

Thay $dV$ từ (15.A1.2), ta có:

Số hạng trong ngoặc đầu tiên (hệ số của $dt$) chính là vế trái của Black-Scholes PDE (15.1), vốn bằng $0$ theo giả thiết $V$ thỏa mãn PDE đó. Vì vậy:

Đây là một phương trình vi phân ngẫu nhiên không có số hạng $dt$, nghĩa là $Y_t$ là một martingale. Tính chất martingale nghĩa là, với $t < T$:

Áp dụng điều kiện cuối $V(S,T)=f(S)$, ta có $Y_T = e^{-rT}V(S_T,T) = e^{-rT}f(S_T)$. Tại thời điểm $t$, ta có $Y_t = e^{-rt}V(S_t,t)$:

Nhân cả hai vế với $e^{rt}$:

Điều cần chứng minh.

A3. Lời giải hệ số $a_i, b_i, c_i$ cho Explicit Method

Thay các xấp xỉ đạo hàm (15.13), tính tại nút $(i,j+1)$, vào Black-Scholes PDE (15.1) cũng tính tại nút $(i,j+1)$:

Thay $S_i = i\Delta S$, số hạng thứ hai và thứ ba trở thành:

Thay số hạng thứ hai và thứ ba vào phương trình (15.A3.1) và nhân toàn bộ với $\Delta t$, ta có:

Nhóm lại theo từng ẩn số cột $j+1$, ta có:

Đặt $a_i$, $b_i$, $c_i$ lần lượt là ba hệ số trên, ta thu được công thức (15.14).

• Hệ số $a_i = -\frac{\Delta t}{2}ri + \frac{\Delta t}{2}\sigma^2 i^2 = \frac{\Delta t}{2}\left(\sigma^2 i^2 - ri\right)$
• Hệ số $b_i = 1 - \Delta t\,\sigma^2 i^2 - \Delta t\, r = 1 - \Delta t\left(\sigma^2 i^2 + r\right)$
• Hệ số $c_i = \frac{\Delta t}{2}ri + \frac{\Delta t}{2}\sigma^2 i^2 = \frac{\Delta t}{2}\left(\sigma^2 i^2 + ri\right)$

Điều cần chứng minh.

A4. Lời giải hệ số $\alpha_i, \beta_i, \gamma_i$ cho Implicit Method

Tương tự Phụ lục (A3), thay các xấp xỉ đạo hàm (15.17), vào Black-Scholes PDE (15.1) bây giờ được tính tại nút $(i,j)$:

Thay $S_i = i\Delta S$, số hạng thứ hai và thứ ba trở thành:

Nhân (15.A4.2) với $\Delta t$ và giải $V_{i,j+1}$ (số hạng duy nhất không thuộc cột $j$):

Nhóm lại theo từng ẩn số của cột $j$:

Đặt $\alpha_i$, $\beta_i$, $\gamma_i$ lần lượt là ba hệ số trên, ta thu được công thức (15.18).

• Hệ số $\alpha_i = \frac{\Delta t}{2} ri - \frac{\Delta t}{2}\sigma^2 i^2 = -\Delta t \frac{1}{2}\left(\sigma^2 i^2 - ri\right)$
• Hệ số $\beta_i = 1 + \Delta t\sigma^2 i^2 + \Delta t r = 1 + \Delta t \left(\sigma^2 i^2 + r\right)$
• Hệ số của $\gamma_i = -\frac{\Delta t}{2} ri - \frac{\Delta t}{2}\sigma^2 i^2 = -\Delta t \frac{1}{2}\left(\sigma^2 i^2 + ri\right)$

Điều cần chứng minh.

A5. Lời giải hệ số cho Crank-Nicolson Method, độ chính xác bậc hai và tính ổn định

Đặt $\mathcal{L}V_{i,k}$ là phần không gian (ba số hạng cuối) của Black-Scholes PDE (15.1), khi đó PDE trở thành $\frac{\partial V}{\partial t} + \mathcal{L}V = 0$. Dùng Central Difference tại cột $k$, và thay $S_i = i\Delta S$, ta có:

Nhóm lại theo từng giá trị $V_{i-1,k}, V_{i,k}, V_{i+1,k}$, ta nhận ra ngay các hệ số của Phụ lục (A3):

Công thức trên đúng cho mọi cột $k$, nên cũng đúng cho cột $j$ và cột $j+1$.

Theo định nghĩa (15.20), Crank-Nicolson giữ nguyên đạo hàm thời gian $\frac{V_{i,j+1}-V_{i,j}}{\Delta t}$ và lấy trung bình cộng đạo hàm theo không gian ở hai cột:

Nhân với $\Delta t$ và thay (15.A5.1) cho cả hai cột $j$ và $j+1$:

Chuyển các ẩn số cột $j$ sang vế trái, các giá trị đã biết cột $j+1$ sang vế phải:

Điều cần chứng minh, công thức (15.21).