利用MATLAB实现以下功能，从外部导入10个合约数据：第一列为内在价值（标的资产价格-执行价），第二列为卖价，第三列为最新价（看作市场价），第四列为执行价，10个合约的到期时间T均为21/365。根据卖价使用Metropolis-Hastings算法对波动率进行反演，要求建议分步为高斯分布，并根据已知的标的资产价格、执行价和到期时间，将反演得到的波动率带回B-S公式，得到期权的理论价格，将理论价格与最新价和卖价进行比较。并计算误差

以下是MATLAB代码实现：

% 导入数据
data = xlsread('contracts.xlsx');
S = data(:,2); % 标的资产价格
K = data(:,4); % 执行价
r = 0.03; % 无风险利率
T = 21/365; % 到期时间
N = length(S); % 合约个数

% 定义Metropolis-Hastings算法参数
niter = 10000; % 迭代次数
sigma0 = 0.2; % 初始波动率
step = 0.02; % 步长
acc = 0; % 接受率

% 定义高斯分布函数
f = @(x) exp(-0.5*x.^2) / sqrt(2*pi);

% Metropolis-Hastings算法反演波动率
sigma = zeros(N,1);
for i = 1:N
    y = log(S(i)/K(i)) + r*T; % 内在价值转化为对数价格
    p = @(s) bs_price(s, K(i), r, T, 0) - data(i,2); % 定义误差函数
    for j = 1:niter
        sigma_new = sigma(i) + step*randn;
        if sigma_new > 0 % 波动率必须为正数
            alpha = min(1, f(p(sigma_new))/f(p(sigma)));
            if rand < alpha % 接受新样本
                sigma(i) = sigma_new;
                acc = acc + 1;
            end
        end
    end
end

% 计算误差
price_bs = zeros(N,1);
for i = 1:N
    price_bs(i) = bs_price(S(i), K(i), r, T, sigma(i));
end
error = abs(price_bs - data(:,3)) ./ data(:,3);

% 显示结果
disp(['Metropolis-Hastings算法接受率为 ', num2str(acc/(N*niter))]);
disp(['平均误差为 ', num2str(mean(error))]);

% 定义B-S模型计算期权价格的函数
function price = bs_price(S, K, r, T, sigma)
    d1 = (log(S/K) + (r + 0.5*sigma^2)*T) / (sigma*sqrt(T));
    d2 = d1 - sigma*sqrt(T);
    price = S*normcdf(d1) - K*exp(-r*T)*normcdf(d2);
end

其中，`xlsread`函数用于从Excel文件中读取数据，`bs_price`函数用于计算B-S模型的期权价格。在程序中，我们首先导入数据，并定义了Metropolis-Hastings算法的参数和高斯分布函数。然后，我们使用循环反演每个合约的波动率，并计算误差。最后，我们显示了反演的结果和误差。