#COPYRIGHT NOTICE:  Copyright (c) 1997-1998 by Sasha Cyganowski.
#                                 2003-2004 by Lars Gruene
#Permission is  
#granted to anyone to use, modify, or redistribute this software freely, 
#subject to the following restrictions:  1. The author accepts no
#responsibility for any consequences of this software and makes no guarantee
#that this software is free of defects.  2.  The origin of this software must  
#not be misrepresented, either by explicit claim or by omission. 
#3.  This notice and the copyright must be included with all copies or
#modified versions of this software.  4.  This software may not be included or
#redistributed as part of any package to be sold for profit unlesss the author
#has given explicit written permission to do so.    
# Sasha Cyganowski, School of Computing and Mathematics, Deakin University,
#Waurn Ponds, VIC, Australia.  Internet: sash@deakin.edu.au
# Telephone:  +61 3 5227 1382.
# additional routines and updates by Lars Gruene (lars.gruene@uni-bayreuth.de)

with(stats):

Wiener:=proc(t,p)\

local temp1,temp2,temp3,temp4,temp5,temp6,i;
temp6[1]:=0;
for i from 1 to p do
temp6[i+1]:=temp6[i]+stats[random,normald[0,1]](1)*convert(sqrt(t),float);
od;

RETURN(convert(temp6,list));
end:


stochplot:=proc(s,h,p,z)
local i,j,temp1,temp2,Xdata,Ydata;

temp1:=Wiener(h,p);

temp2:=rhs(s[1]);
Ydata[1]:=z;
j:=1;

for i from 1 to p do

Ydata[i+1]:=evalf(subs(Y||j[n]=Ydata[i],Delta[n]=h,DeltaW||j[n]=temp1[i+1]-temp1[i],temp2));
od;
Xdata[1]:=0;
for i from 1 to p do
Xdata[i+1]:=Xdata[i]+h;
od;
Xdata:=convert(Xdata,list);
Ydata:=convert(Ydata,list);
plots[display]({statplots[scatterplot](Xdata,Ydata)},axes=FRAME);
end:


stochplot_time:=proc(s,h,p,z,t0)
local i,j,temp1,temp2,Xdata,Ydata;

temp1:=Wiener(h,p);

temp2:=rhs(s[1]);
Ydata[1]:=z;
j:=1;

for i from 1 to p do

Ydata[i+1]:=evalf(subs(Y||j[n]=Ydata[i],Delta[n]=h,DeltaW||j[n]=temp1[i+1]-temp1[i],t=t0+(p-1)*h,temp2));
od;
Xdata[1]:=t0;
for i from 1 to p do
Xdata[i+1]:=Xdata[i]+h;
od;
Xdata:=convert(Xdata,list);
Ydata:=convert(Ydata,list);
plots[display]({statplots[scatterplot](Xdata,Ydata)},axes=FRAME);
end: