Plotting confidence intervals in Matlab

Question

I am having trouble plotting confidence intervals through errorbar function inside Matlab. I have written the following code below

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

in order to match the picture displayed below (for some reason, plot(M,ch08(10,9,0.06,0.1,1),'--') doesn't display anything, but I am ignoring this cosmetic problem).

In the above Matlab code, a confidence interval is calculated by

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

My current implementation almost matches the picture above.

I don't know how to plot confidence intervals inside Matlab. I looked on Google and found that the recommended method is through the errorbar function.

I think that I can add

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

inside the errorbar function in order to plot the vertical confidence interval lines displayed inside the first picture. Is this something that can be accomplished by adjusting

errorbar(aM,conf(1),conf(2))

to somehow track the changes from conf = [aM - 1.96*bM/sqrt(M), aM + 1.96*bM/sqrt(M)]?

I am also referencing a second script in my Matlab code,

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

Answer 1

• errorbar is usually used for this purpose but as you also figured out, line (or plot) can also be used to draw confidence interval lines. I'll focus on using errorbar.
• plot(M,ch08(10,9,0.06,0.1,1),'--.k') doesn't draw dashed lines since ch08(10,9,0.06,0.1,1) is just a value. For a line to be drawn, there needs to be multiple y values (equal to the number of x values) which can be same (just like your case). Otherwise it would be just points (which is happening in your code).

I have incorporated the above and some other optimisations in your code below:

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

Result: