function bpath3 ( )

%*****************************************************************************80
%
%% bpath3() displays the average of 1000 Brownian paths.
%
%  Discussion:
%
%    This routine computes 1000 simulations of discretized Brownian 
%    motion W(t) over the time interval [0,1] using 500 time steps.
%
%    Actually, we are interested in a function u(W(t)):
%
%      u(W(t)) = exp ( t + W(t)/2 )
%
%    The routine plots 5 of the simulations, as well as the average
%    of all the simulations.  
%
%    The plot of the average should be quite smooth.  Its expected
%    value is exp ( 9 * t / 8 ), and we compute the 'error', that is,
%    the difference between the averaged value and this expected
%    value.  This 'error' should decrease as the number of simulation
%    is increased.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 September 2006
%
%  Author:
%
%    Desmond Higham
%
%  Reference:
%
%    Desmond Higham,
%    An Algorithmic Introduction to Numerical Simulation of
%    Stochastic Differential Equations,
%    SIAM Review,
%    Volume 43, Number 3, September 2001, pages 525-546.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'bpath3():\n' );
  fprintf ( 1, '  Average 1000 Brownian path simulations.\n' );

% rng ( 'shuffle' );

  T = 1.0;
  N = 500;
  dt = T / N;
  t = [ dt: dt: 1 ];
%
%  We will be computing M paths simultaneously.
%
  M = 1000;
  dw = sqrt ( dt ) * randn ( M, N );
  w0 = 0.0;
  w = cumsum ( dw, 2 );

  U = exp ( repmat ( t, [M 1] ) + 0.5 * w );

  Umean = mean ( U );
%
%  Plot the mean of the M paths.
%
  plot ( [0,t], [1,Umean], 'b-' );
  hold ( 'on' );
%
%  Overlay the plots of the first 5 paths.
%
  plot ( [0,t], [ones(5,1), U(1:5,:)], 'r--' );
  grid ( 'on' );
  hold ( 'off' );

  xlabel ( 't', 'FontSize', 16 );
  ylabel ( 'U(t)', 'FontSize', 16, 'Rotation', 0, 'HorizontalAlignment', 'right' );
  legend ( 'mean of 1000 paths', '5 individual paths', 2 );
  title ( '1000 paths computed by bpath3()' );

  filename = 'bpath3.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as %s\n', filename );
%
%  Sample error:
%
  averr = norm ( ( Umean - exp(9*t/8) ), 'inf' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum error in averaged data is %f\n', averr );

  return
end
