在 Matlab 中绘制置信区间

Plotting confidence intervals in Matlab

我无法通过 Matlab 中的 errorbar 函数绘制置信区间。我在下面写了下面的代码

clear all;
close all;
%CH15      Program for Chapter 15
%
%   Monte Carlo for a European call

randn('state',100)

%%%%%%%%%%%%%%%%% Problem and method parameters %%%%%%%%%%%%%%%
S = 10; E = 9; sigma = 0.1; r = 0.06; T = 1; 
Dt = 1e-3; N = T/Dt; M = 1e4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hold on;
for M = [2^5,2^6,2^7,2^8,2^9,2^10,2^11,2^12,2^13,2^14,2^15,2^16,2^17]
    V = zeros(M,1);
    for i = 1:M
        Sfinal = S*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*randn);
        V(i) = exp(-r*T)*max(Sfinal-E,0);
    end
    aM = mean(V); bM = std(V);
    conf = [aM - 1.96*bM/sqrt(M), aM + 1.96*bM/sqrt(M)]
    %xlabel('Samples') % x-axis label
    title('Monte Carlo Approximations')
    ylabel('Option value approximation') % y-axis label
    set(gca, 'YScale', 'log')
    set(gca, 'XScale', 'log')
    yticks([10^0.1 10^0.2 10^0.3])
    axis([10^1 10^6 10^0.1 10^0.3])
    yticklabels({'10^{0.1}','10^{0.2}','10^{0.3}'})
    plot(M,aM,'x')
    plot(M,ch08(10,9,0.06,0.1,1),'--.k')
    err = ones*(size(conf));
    errorbar(aM,conf(1),conf(2))
end

为了匹配下面显示的图片(由于某些原因,plot(M,ch08(10,9,0.06,0.1,1),'--') 没有显示任何内容,但我忽略了这个外观问题)。

在上面的Matlab代码中,置信区间的计算方式是

conf = [aM - 1.96*bM/sqrt(M), aM + 1.96*bM/sqrt(M)]

我目前的实现与上图几乎一致。

我不知道如何在 Matlab 中绘制置信区间。我看了Google,发现推荐的方法是通过errorbar函数。

我想我可以添加

conf = [aM - 1.96*bM/sqrt(M), aM + 1.96*bM/sqrt(M)]

在 errorbar 函数中绘制第一张图片中显示的垂直置信区间线。这是可以通过调整

来完成的事情吗?
errorbar(aM,conf(1),conf(2))

以某种方式跟踪 conf = [aM - 1.96*bM/sqrt(M), aM + 1.96*bM/sqrt(M)]?

的变化

我还在我的 Matlab 代码中引用了第二个脚本,

function [C, Cdelta, P, Pdelta] = ch08(S,E,r,sigma,tau)
% Program for Chapter 8
% This is a MATLAB function
%
% Input arguments: S = asset price at time t
%                  E = Exercise price
%                  r = interest rate
%                  sigma = volatility
%                  tau = time to expiry (T-t) 
%
% Output arguments: C = call value, Cdelta = delta value of call 
%                   P = Put value, Pdelta = delta value of put
%
%   function [C, Cdelta, P, Pdelta] = ch08(S,E,r,sigma,tau)

if tau > 0
   d1 = (log(S/E) + (r + 0.5*sigma^2)*(tau))/(sigma*sqrt(tau));
   d2 = d1 - sigma*sqrt(tau);
   N1 = 0.5*(1+erf(d1/sqrt(2)));
   N2 = 0.5*(1+erf(d2/sqrt(2)));
   C = S*N1-E*exp(-r*(tau))*N2;
   Cdelta = N1;
   P = C + E*exp(-r*tau) - S;
   Pdelta = Cdelta - 1;
else 
   C = max(S-E,0);
   Cdelta =  0.5*(sign(S-E) + 1);
   P = max(E-S,0);
   Pdelta = Cdelta - 1;
end

一种解决方案是通过添加

line([M M], conf);

绘制一条垂直直线

clear all;
close all;
%CH15      Program for Chapter 15
%
%   Monte Carlo for a European call

randn('state',100)

%%%%%%%%%%%%%%%%% Problem and method parameters %%%%%%%%%%%%%%%
S = 10; E = 9; sigma = 0.1; r = 0.06; T = 1; 
Dt = 1e-3; N = T/Dt; M = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hold on;
for M = [2^5,2^6,2^7,2^8,2^9,2^10,2^11,2^12,2^13,2^14,2^15,2^16,2^17]
    V = zeros(M,1);
    for i = 1:M
        Sfinal = S*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*randn);
        V(i) = exp(-r*T)*max(Sfinal-E,0);
    end
    aM = mean(V); bM = std(V);
    conf = [aM - 1.96*bM/sqrt(M), aM + 1.96*bM/sqrt(M)]
    title('Monte Carlo Approximations')
    xlabel('Num samples') % x-axis label
    ylabel('Option value approximation') % y-axis label
    set(gca, 'YScale', 'log')
    set(gca, 'XScale', 'log')
    yticks([10^0.1 10^0.2 10^0.3])
    axis([10^1 10^6 10^0.1 10^0.3])
    yticklabels({'10^{0.1}','10^{0.2}','10^{0.3}'})
    plot(M,aM,'x')
    plot(M,ch08(10,9,0.06,0.1,1),'--.k')
    line([M M], conf);
end

在针对

进行测试时,Matlab 中似乎存在错误
plot(M,ch08(10,9,0.06,0.1,1),'--')

不绘制任何数据。

  • errorbar is usually used for this purpose but as you also figured out, line (or plot) 也可用于绘制置信区间线。我将专注于使用 errorbar.
  • plot(M,ch08(10,9,0.06,0.1,1),'--.k') 不绘制虚线,因为 ch08(10,9,0.06,0.1,1) 只是一个值。要绘制一条线,需要有多个 y 值(等于 x 值的数量),这些值可以相同(就像您的情况一样)。否则它只是点数(这在您的代码中发生)。

我已将上述内容和其他一些优化合并到您的以下代码中:

randn('state', 100);
S=10;   E=9;   sigma=0.1;   r=0.06;   T=1;   Dt=1e-3;   N=T/Dt; 
M = [2^5,2^6,2^7,2^8,2^9,2^10,2^11,2^12,2^13,2^14,2^15,2^16,2^17];
hold on;
for k=1:numel(M)
    %No need of loop here. Generate all random values in one go 
    Sfinal = S*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*randn(M(k),1));
    V = exp(-r*T)*max(Sfinal-E,0);
    aM = mean(V);   bM = std(V);

    plot(M(k),aM,'x');
    errorbar(M(k), aM, 1.96*bM/sqrt(M(k)));    
end
chvar = repmat(ch08(10,9,0.06,0.1,1),1,numel(M));  %<----Notice this
plot(M, chvar,'--.k'); 

%Other figure cosmetics
%These commands shouldn't be inside the loop to avoid unnecessary computations
title('Monte Carlo Approximations');
xlabel('Samples'); % x-axis label
ylabel('Option value approximation'); % y-axis label
set(gca,'XScale', 'log','YScale', 'log');
axis([10^1 10^6 10^0.1 10^0.3]);
set(gca,'YTick',[10^0.1 10^0.2 10^0.3]);
set(gca,'YTickLabel',{'10^{0.1}','10^{0.2}','10^{0.3}'});
%Nothing is wrong with using yticks and yticklabels function but those require >=R2016b

结果: