%
% calculation of multi-dimensional Brownian Bridge
%
% function dW = bb(Z,T)
% function dW = bb(Z,T,C)
%
% T       -- time interval
% Z(:,:)  -- input vectors of unit Normal variables
% C(:,:)  -- square PCA matrix for multi-dimensional Brownian motion (optional)
% dW(:,:) -- output vectors of Brownian path increments
%
% first dimension of Z/dW corresponds to Brownian dimension * #timesteps
% second dimension of Z/dW corresponds to #paths
%

function W = bb(W,T,C)

if nargin<2 | nargin>3
  error('wrong number of arguments')

%
% one-dimensional Brownian motion
%

elseif nargin==2
  K = 1;
  N = size(W,1);

%
% for multi-dimensional Brownian motion, reorder input
% to get contiguous time
%

elseif nargin==3
  K = size(C,1);
  if size(C,1)~=K
    error('error: matrix C not square')
  end
  N = size(W,1);
  if mod(N,K)~=0
    error('error: first dimension of W not divisible by first dimension of C')
  end
  N = N/K;

  p = reshape((reshape((1:N*K)',K,N))',N*K,[]);
  W = reshape(W(p,:),N,[]);
end

%
% perform Brownian Bridge interpolation
%

M = round(log(N)/log(2));
if N~=2^M
  error('error: number of timesteps not a power of 2')
end

for m = 1:M
  W([1:2:2^m 2:2:2^m],:) = ...
   [ 0.5*W(1:2^(m-1),:)+W(2^(m-1)+1:2^m,:)/sqrt(2)^(m+1) ; ...
     0.5*W(1:2^(m-1),:)-W(2^(m-1)+1:2^m,:)/sqrt(2)^(m+1) ];
end

%
% for multi-dimensional Brownian motion, reorder output 
% to get contiguous multiple dimensions and apply PCA
%

if nargin==3
  W = reshape(W,N*K,[]);
  p = reshape((reshape((1:N*K)',N,K))',N*K,[]);
  W = reshape(W(p,:),K,[]);
  W = reshape(C*W,N*K,[]);
end

%
% finally, adjust for non-unit time interval
%

W = sqrt(T)*W;

return

