BSDE（後向隨機微分方程）的歐拉算法實現。

BSDE在金融數學、隨機控制等領域有重要應用。

1. BSDE基本理論

BSDE的一般形式

BSDE通常表示為：

Y_t = ξ + ∫_t^T f(s, Y_s, Z_s) ds - ∫_t^T Z_s dW_s

其中：

  • Y_t 是狀態過程
  • Z_t 是控制過程
  • W_t 是布朗運動
  • ξ 是終端條件
  • f 是生成器函數

2. 歐拉離散化格式

時間離散化

將時間區間 [0, T] 離散為：

0 = t_0 < t_1 < ... < t_N = T

步長：Δt = T/N

歐拉格式

前向-後向歐拉格式：

Y_N = ξ
Z_N = 0

對於 i = N-1, ..., 0:
    Z_i = (1/Δt) E[Y_{i+1} ΔW_i | F_{t_i}]
    Y_i = E[Y_{i+1} | F_{t_i}] + Δt f(t_i, Y_i, Z_i)

3. MATLAB實現

基本BSDE求解器

classdef BSDESolver < handle
    % BSDE求解器 - 歐拉方法實現
    
    properties
        T          % 終端時間
        N          % 時間步數
        dt         % 時間步長
        M          % 蒙特卡洛路徑數
        d          % 布朗運動維度
        t_grid     % 時間網格
    end
    
    methods
        function obj = BSDESolver(T, N, M, d)
            % 構造函數
            obj.T = T;
            obj.N = N;
            obj.M = M;
            obj.d = d;
            obj.dt = T / N;
            obj.t_grid = linspace(0, T, N+1);
        end
        
        function [Y, Z, paths] = solve(obj, f, phi, X0)
            % 求解BSDE
            % 輸入：
            %   f: 生成器函數 f(t, x, y, z)
            %   phi: 終端條件函數 phi(x)
            %   X0: 初始狀態
            
            % 初始化
            Y = zeros(obj.M, obj.N+1);
            Z = zeros(obj.M, obj.d, obj.N+1);
            X = zeros(obj.M, obj.N+1);
            W = zeros(obj.M, obj.d, obj.N+1);
            
            % 生成布朗運動路徑
            X(:,1) = X0;
            for i = 1:obj.N
                dW = sqrt(obj.dt) * randn(obj.M, obj.d);
                W(:,:,i+1) = W(:,:,i) + dW;
                % 這裏可以添加前向過程的具體演化
                X(:,i+1) = X(:,i) + dW; % 簡單布朗運動示例
            end
            
            % 設置終端條件
            Y(:,end) = phi(X(:,end));
            Z(:,:,end) = 0;
            
            % 後向迭代
            for i = obj.N:-1:1
                % 計算條件期望
                [E_Y, E_YZ] = obj.conditional_expectation(Y(:,i+1), W(:,:,i+1));
                
                % 更新Z過程
                dW = W(:,:,i+1) - W(:,:,i);
                for j = 1:obj.d
                    Z(:,j,i) = E_YZ(:,j) / obj.dt;
                end
                
                % 更新Y過程（隱式歐拉）
                Y_old = E_Y;
                max_iter = 10;
                tol = 1e-6;
                
                for iter = 1:max_iter
                    Y_new = E_Y + obj.dt * f(obj.t_grid(i), X(:,i), Y_old, Z(:,:,i));
                    if max(abs(Y_new - Y_old)) < tol
                        break;
                    end
                    Y_old = Y_new;
                end
                Y(:,i) = Y_new;
            end
            
            paths.X = X;
            paths.W = W;
        end
        
        function [E_Y, E_YZ] = conditional_expectation(obj, Y, W)
            % 計算條件期望
            % 使用最小二乘蒙特卡洛方法
            
            E_Y = zeros(obj.M, 1);
            E_YZ = zeros(obj.M, obj.d);
            
            % 使用當前狀態作為基函數
            basis = [ones(obj.M,1), W];
            
            % 迴歸計算E[Y|F_t]
            coeff_Y = (basis' * basis) \ (basis' * Y);
            E_Y = basis * coeff_Y;
            
            % 迴歸計算E[YΔW|F_t]
            for j = 1:obj.d
                Y_dW = Y .* W(:,j);
                coeff_YZ = (basis' * basis) \ (basis' * Y_dW);
                E_YZ(:,j) = basis * coeff_YZ;
            end
        end
    end
end

具體的BSDE示例實現

function bsde_euler_demo()
    % BSDE歐拉算法演示
    
    %% 參數設置
    T = 1.0;        % 終端時間
    N = 100;        % 時間步數
    M = 10000;      % 蒙特卡洛路徑數
    d = 1;          % 布朗運動維度
    X0 = 1.0;       % 初始狀態
    
    %% 創建求解器
    solver = BSDESolver(T, N, M, d);
    
    %% 示例1: 線性BSDE
    fprintf('求解線性BSDE...\n');
    [Y_lin, Z_lin, paths_lin] = solve_linear_bsde(solver, X0);
    
    %% 示例2: 非線性BSDE
    fprintf('求解非線性BSDE...\n');
    [Y_nonlin, Z_nonlin, paths_nonlin] = solve_nonlinear_bsde(solver, X0);
    
    %% 可視化結果
    plot_bsde_results(Y_lin, Z_lin, Y_nonlin, Z_nonlin, paths_lin, paths_nonlin);
end

function [Y, Z, paths] = solve_linear_bsde(solver, X0)
    % 線性BSDE示例
    % dY_t = (aY_t + bZ_t) dt + Z_t dW_t
    % Y_T = φ(X_T)
    
    a = 0.1;
    b = 0.2;
    
    % 生成器函數
    f = @(t, x, y, z) a * y + b * z;
    
    % 終端條件
    phi = @(x) x.^2;  % Y_T = X_T^2
    
    [Y, Z, paths] = solver.solve(f, phi, X0);
end

function [Y, Z, paths] = solve_nonlinear_bsde(solver, X0)
    % 非線性BSDE示例
    % dY_t = (Y_t^2 + Z_t^2) dt + Z_t dW_t
    % Y_T = φ(X_T)
    
    % 非線性生成器函數
    f = @(t, x, y, z) y.^2 + z.^2;
    
    % 終端條件
    phi = @(x) sin(x);  % Y_T = sin(X_T)
    
    [Y, Z, paths] = solver.solve(f, phi, X0);
end

4. 改進的歐拉算法

帶有投影步的歐拉方法

classdef ProjectedBSDESolver < BSDESolver
    % 帶投影的BSDE求解器
    
    methods
        function [Y, Z, paths] = solve_with_projection(obj, f, phi, X0, projection)
            % 帶投影的BSDE求解
            % projection: 投影函數
            
            [Y, Z, paths] = solve@BSDESolver(obj, f, phi, X0);
            
            % 應用投影
            for i = 1:obj.N+1
                Y(:,i) = projection(Y(:,i));
                for j = 1:obj.d
                    Z(:,j,i) = projection(Z(:,j,i));
                end
            end
        end
    end
end

高階歐拉方法

classdef HighOrderBSDESolver < BSDESolver
    % 高階BSDE求解器
    
    methods
        function [Y, Z, paths] = solve_high_order(obj, f, phi, X0)
            % 高階歐拉方法
            
            % 初始化
            Y = zeros(obj.M, obj.N+1);
            Z = zeros(obj.M, obj.d, obj.N+1);
            X = zeros(obj.M, obj.N+1);
            
            X(:,1) = X0;
            
            % 生成路徑
            for i = 1:obj.N
                dW = sqrt(obj.dt) * randn(obj.M, obj.d);
                X(:,i+1) = X(:,i) + dW;
            end
            
            % 終端條件
            Y(:,end) = phi(X(:,end));
            
            % 高階後向迭代
            for i = obj.N:-1:1
                % 預測步
                [E_Y_pred, E_YZ_pred] = obj.high_order_expectation(Y(:,i+1), X(:,i));
                
                % 校正步
                Y_pred = E_Y_pred + obj.dt * f(obj.t_grid(i), X(:,i), E_Y_pred, E_YZ_pred/obj.dt);
                [E_Y_corr, E_YZ_corr] = obj.high_order_expectation(Y_pred, X(:,i));
                
                % 更新
                for j = 1:obj.d
                    Z(:,j,i) = E_YZ_corr(:,j) / obj.dt;
                end
                Y(:,i) = 0.5 * (E_Y_pred + E_Y_corr) + obj.dt * f(obj.t_grid(i), X(:,i), Y_pred, Z(:,:,i));
            end
            
            paths.X = X;
        end
        
        function [E_Y, E_YZ] = high_order_expectation(obj, Y, X)
            % 高階條件期望計算
            % 使用多項式基函數
            
            E_Y = zeros(obj.M, 1);
            E_YZ = zeros(obj.M, obj.d);
            
            % 高階基函數
            basis = [ones(obj.M,1), X, X.^2, X.^3];
            
            % 迴歸計算條件期望
            coeff_Y = (basis' * basis) \ (basis' * Y);
            E_Y = basis * coeff_Y;
            
            % 為每個布朗運動維度計算
            for j = 1:obj.d
                Y_basis = Y .* basis(:,2);  % 使用X作為權重
                coeff_YZ = (basis' * basis) \ (basis' * Y_basis);
                E_YZ(:,j) = basis * coeff_YZ;
            end
        end
    end
end

5. 金融應用示例

歐式期權定價的BSDE

function black_scholes_bsde_demo()
    % Black-Scholes模型的BSDE求解
    
    T = 1.0;
    N = 100;
    M = 5000;
    d = 1;
    
    % 市場參數
    S0 = 100;       % 初始股價
    K = 100;        % 行權價
    r = 0.05;       % 無風險利率
    sigma = 0.2;    % 波動率
    
    solver = BSDESolver(T, N, M, d);
    
    % Black-Scholes生成器函數
    f = @(t, s, y, z) -r * y + (r / (sigma * s)) * z;
    
    % 看漲期權終端條件
    phi = @(s) max(s - K, 0);
    
    [Y, Z, paths] = solver.solve(f, phi, S0);
    
    % 與解析解比較
    [call_analytical, ~] = blsprice(S0, K, r, T, sigma);
    call_bsde = mean(Y(:,1));
    
    fprintf('BSDE方法期權價格: %.4f\n', call_bsde);
    fprintf('解析解期權價格: %.4f\n', call_analytical);
    fprintf('相對誤差: %.4f%%\n', abs(call_bsde - call_analytical)/call_analytical*100);
    
    plot_option_results(Y, Z, paths, call_analytical);
end

function plot_option_results(Y, Z, paths, analytical_price)
    % 繪製期權定價結果
    
    figure('Position', [100, 100, 1200, 800]);
    
    % 價格過程
    subplot(2,2,1);
    plot(paths.t_grid, Y(1:10,:)');
    hold on;
    plot(paths.t_grid, mean(Y), 'k-', 'LineWidth', 3);
    title('期權價格過程 Y_t');
    xlabel('時間');
    ylabel('價格');
    legend('樣本路徑', '平均值', 'Location', 'best');
    
    % Delta對衝
    subplot(2,2,2);
    plot(paths.t_grid, squeeze(Z(1:10,1,:))');
    hold on;
    plot(paths.t_grid, mean(squeeze(Z(:,1,:))), 'k-', 'LineWidth', 3);
    title('Delta對衝係數 Z_t');
    xlabel('時間');
    ylabel('Delta');
    
    % 價格分佈
    subplot(2,2,3);
    histogram(Y(:,1), 50, 'Normalization', 'probability');
    hold on;
    line([analytical_price, analytical_price], ylim, 'Color', 'red', 'LineWidth', 2);
    title('初始價格分佈');
    xlabel('價格');
    ylabel('概率密度');
    legend('BSDE價格', '解析價格');
    
    % 收斂性分析
    subplot(2,2,4);
    time_points = 1:size(Y,2);
    mean_prices = mean(Y);
    std_prices = std(Y);
    plot(paths.t_grid, mean_prices, 'b-', 'LineWidth', 2);
    hold on;
    fill([paths.t_grid, fliplr(paths.t_grid)], ...
         [mean_prices + std_prices, fliplr(mean_prices - std_prices)], ...
         'b', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
    title('價格收斂過程');
    xlabel('時間');
    ylabel('價格');
end

6. 收斂性分析

function convergence_analysis()
    % BSDE歐拉方法的收斂性分析
    
    T = 1.0;
    M = 10000;
    d = 1;
    X0 = 1.0;
    
    % 不同時間步數
    N_values = [25, 50, 100, 200, 400];
    errors_Y = zeros(size(N_values));
    errors_Z = zeros(size(N_values));
    
    % 參考解（使用最精細的網格）
    N_ref = 1000;
    solver_ref = BSDESolver(T, N_ref, M, d);
    f_ref = @(t, x, y, z) 0.1 * y + 0.2 * z;
    phi_ref = @(x) x.^2;
    [Y_ref, Z_ref, ~] = solver_ref.solve(f_ref, phi_ref, X0);
    Y0_ref = mean(Y_ref(:,1));
    Z0_ref = mean(Z_ref(:,1,1));
    
    for i = 1:length(N_values)
        N = N_values(i);
        solver = BSDESolver(T, N, M, d);
        [Y, Z, ~] = solver.solve(f_ref, phi_ref, X0);
        
        Y0 = mean(Y(:,1));
        Z0 = mean(Z(:,1,1));
        
        errors_Y(i) = abs(Y0 - Y0_ref);
        errors_Z(i) = abs(Z0 - Z0_ref);
        
        fprintf('N = %d, Y0誤差 = %.6f, Z0誤差 = %.6f\n', ...
                N, errors_Y(i), errors_Z(i));
    end
    
    % 繪製收斂圖
    figure;
    loglog(N_values, errors_Y, 'o-', 'LineWidth', 2, 'DisplayName', 'Y誤差');
    hold on;
    loglog(N_values, errors_Z, 's-', 'LineWidth', 2, 'DisplayName', 'Z誤差');
    loglog(N_values, 1./N_values, '--', 'DisplayName', 'O(1/N)');
    loglog(N_values, 1./sqrt(N_values), '--', 'DisplayName', 'O(1/√N)');
    
    xlabel('時間步數 N');
    ylabel('誤差');
    title('BSDE歐拉方法收斂性');
    legend('Location', 'best');
    grid on;
end

7. 使用示例

% 基本使用示例
clear; clc; close all;

% 運行演示
bsde_euler_demo();

% 金融應用示例
black_scholes_bsde_demo();

% 收斂性分析
convergence_analysis();

% 自定義BSDE求解
T = 1.0; N = 100; M = 5000; d = 1; X0 = 1.0;
solver = BSDESolver(T, N, M, d);

% 定義自己的生成器和終端條件
my_f = @(t, x, y, z) -0.5 * y.^2 + z;  % 非線性生成器
my_phi = @(x) exp(-0.1 * x.^2);        % 終端條件

[Y, Z, paths] = solver.solve(my_f, my_phi, X0);

fprintf('初始值 Y0 = %.6f\n', mean(Y(:,1)));

參考代碼 解BSDE方程的歐拉算法 www.youwenfan.com/contentcnm/80343.html

算法特點總結

方法

收斂階

計算複雜度

適用場景

基本歐拉

O(Δt)

O(M×N)

一般BSDE

投影歐拉

O(Δt)

O(M×N)

帶約束問題

高階歐拉

O(Δt²)

O(M×N)

高精度需求

BSDE歐拉算法實現提供了：

  • 基本歐拉方法
  • 改進的高階方法
  • 金融應用示例
  • 收斂性分析
  • 可視化工具