function Monte-Carlo1( )
  N=50000;
  NRepl=50000;
  r0=0.05;
  mu=0.065;
  sigma=0.02;
  k=0.2;
  P=100000;
  b=0.03;
  VO(N,P,r0,b,k,mu,sigma,NRepl);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function vo=VO(N,P,r0,b,k,mu,sigma,NRepl);
vrpaths=vasic1(r0,k,mu,sigma,N,NRepl);
vo=zeros(NRepl,N);
vo(:,1)=0;
A=P*(r0*(1+r0)^N/((1+r0)^N-1));
for i=1:NRepl
for j=1:N-1
PV(j)=(1+b)*P*(((1+r0)^N-(1+r0)^j)/((1+r0)^N-1));
AV=0;
for t=1:N-j
av=A/(1+vrpaths(i,j+1))^t;
AV=AV+av;
end
MV=AV-PV(j);
dis(1)=1;
for g=1:j
dis(g+1)=dis(g)/(1+vrpaths(i,j));
end
if MV>=0
vo(i,j+1)=MV*dis(j+1);
else
vo(i,j+1)=0;
end
end
end
for i=1:NRepl;
vod(i)=max(vo(i,:));
end
[h,stats]=cdfplot(vod) % statistic statement

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function vrpaths=vasic1(r0,k,mu,sigma,N,NRepl)
vrpaths=zeros(NRepl,N);
vrpaths(:,1)=r0;
rand('seed',0)
for i=1:NRepl
for j=1:N-1
vrpaths(i,j+1)=vrpaths(i,j)+k*(mu-vrpaths(i,j))+sigma*randn;
end
end