clear;
clear
N=200;
dx=1/N;
a=2/pi;
b=2/pi^2;
x=[0:dx:1];
y=[0:dx:1];
M=80;
z=zeros(N+1);
ksi=randn(N+1);
for m=1:M
    for k=1:M
        %z=z+a*sin(pi*k*x')*sin(pi*m*y)*ksi(m,k)/sqrt(m^2+k^2); % GFF
        %z=z+b*sin(pi*k*x')*sin(pi*m*y)*ksi(m,k)/(m*k); %Tensor of two BB
        %z=z+b*sin(pi*(k-0.5)*x')*sin(pi*m*y)*ksi(m,k)/(m*(k-0.5)); %Kiefer: Tensor of BM and BB
        z=z+b*sin(pi*(k-0.5)*x')*sin(pi*(m-0.5)*y)*ksi(m,k)/((m-0.5)*(k-0.5)); %Br. Sheet: Tensor of BM and BM
    end
end
mesh(x,y,z)
hold
%title('GFF on [0,1]x[0,1]; 80^2 terms, 200^2 points'); xlabel('X-axis'); ylabel('Y-axis');
%title('BB-BB on [0,1]x[0,1]; 80^2 terms, 200^2 points'); xlabel('X-axis'); ylabel('Y-axis');
%title('Kiefer Field on [0,1]x[0,1]; 80^2 terms, 200^2 points'); xlabel('BrM'); ylabel('BrBr');
title('Brownian Sheet on [0,1]x[0,1]; 80^2 terms, 200^2 points'); xlabel('X-axis'); ylabel('Y-axis');


% using the argument about crossing a linear barrier, we concluded that E sup_{t>0} W(t)/(1+t) = \sqrt{pi/8} \approx 0.63
% Monte Carlo simulations confirm it. Note that if X(t)=W(t)/(1+t), then dX(t)=-X(t)dt/(1+t) + dW(t)/(1+t). 

clear; dt=0.001; ddt=sqrt(dt); M=200; N=10000; ksi=randn(N,M); for n=1:M X(1,n)=0; for k=1:N-1 X(k+1,n)=X(k,n)+(ddt*ksi(k,n)-dt*X(k,n))/(1+(k+1)*dt); end; end; sum(max(X))/M