% Vasicek process
% dX = a(c - X) dt + b dW 
% compute the probability of the process going over a certain level 

a = 1 ;
b = 4 ; 
c = 8 ; 

T = 5 ;
N = 500 ;
h = T/N ;
t = (0:N)*h ; 

M = 10000 ; 

X0 = 0 ; 
X = X0*ones(M,N+1);
Z = randn(M,N); 

for i=1:N
    X(:,i+1) = X(:,i) + a*(c-X(:,i))*h + b * sqrt(h)*Z(:,i) ;
end

p = sum( max(X,[],2)>11  )/M 
err = sqrt(p*(1-p))/sqrt(M)  