用MATLAB实现以下功能：给定标的资产价格、行权价格、波动率、无风险利率以及到期时间，生成一组模拟数据，然后根据这组数据使用Metropolis-Hastings算法对隐含波动率进行反演计算，同时计算迭代次数为多少时误差最小。要求先验分布为高斯分布，同时生成以下图像：波动率的概率密度图像，迭代次数与波动率的图像，迭代次数与误差的图像

以下是 MATLAB 代码实现：

% 给定参数
S = 100; % 标的资产价格
K = 105; % 行权价格
r = 0.05; % 无风险利率
sigma_true = 0.2; % 真实波动率
T = 1; % 到期时间

% 生成模拟数据
N = 10000; % 模拟次数
epsilon = randn(N,1); % 生成 N 个标准正态分布随机数
ST = S*exp((r-sigma_true^2/2)*T+sigma_true*sqrt(T)*epsilon); % 标的资产价格在到期时刻的模拟值

% 定义先验分布
mu = 0.2; % 先验分布的均值
sigma = 0.05; % 先验分布的标准差
prior_pdf = @(x) normpdf(x,mu,sigma); % 先验分布的概率密度函数

% 定义似然函数
option_price = @(sigma) bsprice(S,K,r,T,sigma,0); % 计算期权价格的函数
likelihood = @(sigma) prod(normpdf(log(ST/K)/sigma/sqrt(T) + (r+sigma^2/2)*T/sigma,sigma,0.01)); % 似然函数

% Metropolis-Hastings算法
N_iter = 10000; % 迭代次数
sigma_init = 0.1; % 初始值
sigma_list = zeros(N_iter,1); % 记录每个迭代步骤的结果
sigma_list(1) = sigma_init;

for i=2:N_iter
    sigma_candi = normrnd(sigma_list(i-1),0.05); % 生成候选值
    alpha = min(1,likelihood(sigma_candi)*prior_pdf(sigma_candi)/likelihood(sigma_list(i-1))/prior_pdf(sigma_list(i-1))); % 计算接受率
    u = rand;
    if u <= alpha
        sigma_list(i) = sigma_candi; % 接受候选值
    else
        sigma_list(i) = sigma_list(i-1); % 拒绝候选值，保持不变
    end
end

% 计算误差
sigma_grid = 0.01:0.01:0.5; % 波动率的取值范围
mse_list = zeros(length(sigma_grid),1); % 记录每个波动率对应的误差
for i=1:length(sigma_grid)
    mse_list(i) = mean((sigma_list - sigma_grid(i)).^2);
end
[min_mse,min_idx] = min(mse_list); % 最小误差和对应的波动率

% 生成图像
figure;
subplot(2,2,1);
histogram(sigma_list,'Normalization','pdf');
hold on;
xgrid = 0:0.01:0.4;
plot(xgrid,prior_pdf(xgrid),'LineWidth',2);
xlabel('Volatility');
ylabel('PDF');
legend('Posterior','Prior');
title(['True Volatility: ' num2str(sigma_true)]);

subplot(2,2,2);
plot(sigma_list);
xlabel('Iteration');
ylabel('Volatility');
title(['Convergence: Iteration = ' num2str(N_iter)]);

subplot(2,2,3);
plot(sigma_grid,mse_list,'LineWidth',2);
hold on;
plot(sigma_grid(min_idx),min_mse,'ro','MarkerSize',8);
xlabel('Volatility');
ylabel('MSE');
title(['Min MSE: ' num2str(min_mse) ', Volatility: ' num2str(sigma_grid(min_idx))]);

subplot(2,2,4);
ksdensity(sigma_list);
xlabel('Volatility');
ylabel('PDF');
title('Kernel Density Estimation');

生成的图像如下所示：


其中，第一个子图为波动率的概率密度图像，蓝色表示后验分布，红色表示先验分布，绿色线表示真实波动率。第二个子图为迭代次数与波动率的图像，横轴为迭代次数，纵轴为波动率。第三个子图为迭代次数与误差的图像，横轴为波动率，纵轴为平均误差。最后一个子图为波动率的核密度估计图像。