* 1.1.0  Apr 23, '97  (Jeroen Weesie/ICS)  STB-41 tt7

set output error
version 2.1
noisily tut_chk randwalk
set display page 23
set more 1
drop _all

#del ;

di _n(2) 
in wh _dup(75) "-" in gr _n
"Tutorial on random walks" _skip(21) "(Albert Verbeek/Jeroen Weesie)"_n(2)
"This tutorial displays some discrete time random walks in 1 and 2"  _n
"dimensions, illustrating the use of " in wh "graph " 
in gr "and the use of random numbers. If" _n
"you run the demo more than once, you will get different graphs each time."
_n(2)
"First a one-dimensional, normal random walk" _n
in wh _dup(75) "-" in gr _n ;

di in wh ". set obs 500 " _sk(39) "/* 500 observations */" ;
push set obs 500 ;
nois set obs 500 ;

di in wh ". gen int t = _n " _sk(44) "/* t = time */ " ;
push gen int t = _n ;
nois gen int t = _n ;

di in wh ". lab var t time = _n"  ;                                            
push lab var t "time = _n" ;
nois lab var t "time = _n" ;

di in wh ". gen sumz = sum(invnorm(uniform())) " _sk(9) 
         "/* to be explained shortly */"; 
push gen sumz = sum(invnorm(uniform())) ;    
nois gen sumz = sum(invnorm(uniform())) ;    

di _n in wh _dup(75) "-" _n in gr
"Explanation:" _n(2)
in wh "uniform() " in gr 
"generates a random number from the uniform distribution on (0,1)." _n
"If it is transformed by the inverse of the cumulative distribution of a" _n
"random variable X, the result will have the same distribution as X. Here," _n
in wh "invnorm() " 
in gr "is the inverse of the cumulative standard normal distribution. " _n
"So " in wh "invnorm(uniform())" 
in gr " is a standard-normally distributed random number." _n
"By taking the " in wh "sum()" 
in gr ", we get a random walk: a sum of independently and " _n
"identically distributed random numbers." _n
in wh _dup(75) "-" _n ;

set more 0; more ; set more 1 ;

di in wh _dup(75) "-" _n in gr
"Next we will plot the random walk against time (=" in wh "_n" in gr ")." _n
in wh _dup(75) "-" _n ;
   
di in wh ". graph sumz t, c(l) s(i) ti(...) yline(0)" ;
push graph sumz t, 
     c(l) s(i) ti("A one-dimensional, normal random walk sum(z)") yline(0) ;
nois graph sumz t, 
     c(l) s(i) ti("A one-dimensional, normal random walk sum(z)") yline(0) ; 

set more 0; more ; set more 1 ;

di _n in wh _dup(75) "-" _n in gr
"Next the sample path of the means" _n 
in wh _dup(75) "-" _n ;

di in wh ". gen meanz = sumz/_n" ;
push gen meanz = sumz/_n ;
nois gen meanz = sumz/_n ;

di in wh ". graph meanz t, c(l) s(i) ti(...) yline(0)" ; 
push graph meanz t, 
     c(l) s(i) ti("The sample path of the means sum(z)/_n") yline(0) ;
nois graph meanz t, 
     c(l) s(i) ti("The sample path of the means sum(z)/_n") yline(0) ;

set more 0; more ; set more 1 ;

di _n in wh _dup(75) "-" _n in gr
"The means converge to the true mean, which is zero. The standard error " _n
"is 1/sqrt(n). So if we multiply the means by " in wh "sqrt(_n)" 
in gr ", we produce a random " _n 
"walk with expectation 0 and constant variance:" _n
in wh _dup(75) "-" _n ;
   
di in wh ". gen z = meanz*sqrt(_n)" ;
push gen z = meanz*sqrt(_n) ;
nois gen z = meanz*sqrt(_n) ;

di in wh ". graph z t, c(l) s(i) ti(...) yline(0)" ;
push graph z t, c(l) s(i) ti("The sample path of meanz*sqrt(_n)") yline(0); 
nois graph z t, c(l) s(i) ti("The sample path of meanz*sqrt(_n)") yline(0); 

set more 0; more ; set more 1 ;

di _n in wh _dup(75) "-" _n in gr
"Now the same graphs, overlayed with a Cauchy random walk. Recall that the " _n
"Cauchy distribution has no expectation, and has  infinite variance. " _n
"Moreover the mean of n independently, Cauchy  distributed variates is " _n
"Cauchy distributed again. Cauchy distributions have very Xheavy tailsX:" _n 
"Cauchy samples often have a few extreme outliers, dominating the whole " _n 
"picture. " _n(2)
"We can generate a Cauchy random variate as the quotient of two independent" _n
"normal variates. " _n
in wh _dup(75) "-" _n ;

di in wh ". gen sumC = sum(invnorm(uniform())/invnorm(uniform()))";
push gen sumC = sum(invnorm(uniform())/invnorm(uniform())) ;
nois gen sumC = sum(invnorm(uniform())/invnorm(uniform())) ;

di in wh ". graph sumz sumC t, c(ll) s(ii) ti(...) yline(0)" ;  
push graph sumz sumC t, c(ll) s(ii) ti("One-dimensional normal and Cauchy random walks") yline(0) ;
nois graph sumz sumC t, c(ll) s(ii) ti("One-dimensional normal and Cauchy random walks") yline(0) ; 

set more 0; more ; set more 1 ;
   
di _n in wh _dup(75) "-" _n in gr
"Next we graph the sample path of the means." _n 
in wh _dup(75) "-" _n ;
   
di in wh ". gen meanC = sumC/_n" ;
push gen meanC = sumC/_n ;
nois gen meanC = sumC/_n ;

di in wh ". graph meanz meanC t, c(ll) s(ii) ti(...) yline(0)" ; 
push graph meanz meanC t, c(ll) s(ii) ti("The sample path of the normal and Cauchy means") yline(0) ;
nois graph meanz meanC t, c(ll) s(ii) ti("The sample path of the normal and Cauchy means") yline(0) ;

set more 0; more ; set more 1 ;

di _n in wh _dup(75) "-" _n in gr
"Explanation: Each extreme outlier dominates the sum for a long time. That" _n
"is, for a long time, the sum is aprroximately constant. There, the path " _n
"of the sample means is approximately proportional to 1/t, displaying part" _n
"of a hyperbola. " _n(2)
"Multiply the mean by " in wh "sqrt(_n)" in gr " to stabilize variance" _n 
in wh _dup(75) "-" _n ;

di in wh ". gen C = meanC*sqrt(_n)" ;
push gen C = meanC*sqrt(_n) ;
nois gen C = meanC*sqrt(_n) ;

di in wh ". graph z C t, c(ll) s(ii) ti(...) yline(0)" ; 
push graph z C t, c(ll) s(ii) ti("The sample path means*sqrt(n)") yline(0) ;
nois graph z C t, c(ll) s(ii) ti("The sample path means*sqrt(n)") yline(0) ;

set more 0; more ; set more 1 ;

di _n in wh _dup(75) "-" _n in gr
"Next a two-dimensional, normal random walk = a Brownian motion, the motion" _n
"of a particle of pollen in water as seen through a microscope. It was " _n
"first observed by Robert Brown in 1827, and explained by Albert Einstein" _n
"in 1905. The motion is caused by the random collisions with the water " _n
"molecules. It is visible evidence of the existence of molecules and of " _n
"their movement. Einstein also described the basic properties of normal " _n
"random walks." _n
in wh _dup(75) "-" _n ;

di in wh ". gen sumz2 = sum(invnorm(uniform()))" ;
push gen sumz2 = sum(invnorm(uniform())) ;
nois gen sumz2 = sum(invnorm(uniform())) ;

di in wh ". graph sumz sumz2, c(l) s(i) ti(...) yline(0)  xline(0)  border b2( )";
push graph sumz sumz2, c(l) s(i) ti("Two dimensional normal random walk, or Brownian motion") yline(0)  xline(0)  border b2(" ");
nois graph sumz sumz2, c(l) s(i) ti("Two dimensional normal random walk, or Brownian motion") yline(0)  xline(0)  border b2(" ");

set more 0; more ; set more 1 ;

di _n in wh _dup(75) "-" _n in gr
"Some interesting properties of two-dimensional, normal random walks:" _n
"  - At time t the expected distance from the origin is sqrt(t)." _n
"  - Every open disk in the plane, however small, has probability one of" _n
"    being visited by the random walk (if it continues for ever). A similar" _n 
"    statement is NOT true for 3 dimensional walks." _n
"  - The normal random walk is rotationally invariant; there are no" _n
"    " _quote "preferred directions." _quote _n(2)
"Here is another normal walk in the plane." _n
in wh _dup(75) "-" _n ;
    
di in wh ". replace sumz  = sum(invnorm(uniform()))";
push replace sumz  = sum(invnorm(uniform()));
nois replace sumz  = sum(invnorm(uniform()));

di in wh ". replace sumz2 = sum(invnorm(uniform()))";
push replace sumz2 = sum(invnorm(uniform()));
nois replace sumz2 = sum(invnorm(uniform()));

di in wh 
". graph sumz sumz2, c(l) s(i) ti(...) yline(0) xline(0) border b2( )";
push graph sumz sumz2, 
     c(l) s(i) ti("Two dimensional normal random walk, or Brownian motion") 
     yline(0) xline(0) border b2(" ") ;
nois graph sumz sumz2, 
     c(l) s(i) ti("Two dimensional normal random walk, or Brownian motion") 
     yline(0) xline(0) border b2(" ") ;

set more 0; more ; set more 1 ;

di _n in wh _dup(75) "-" _n in gr
"Next we will show a Cauchy random walk in the plane. Can you predict" _n
"how it will differ qualitatively from a normal random walk?" _n 
in wh _dup(75) "-" _n ;

di in wh ". gen sumC2   = sum(invnorm(uniform())/invnorm(uniform()))" ;
push gen sumC2   = sum(invnorm(uniform())/invnorm(uniform())) ;
nois gen sumC2   = sum(invnorm(uniform())/invnorm(uniform())) ;

di in wh 
". graph sumC sumC2, c(l) s(i) ti(...) yline(0) xline(0) border b2( )";
push graph sumC sumC2, c(l) s(i) ti("Two dimensional Cauchy random walk") 
     		       yline(0) xline(0) border b2(" ");
nois graph sumC sumC2, c(l) s(i) ti("Two dimensional Cauchy random walk") 
     		       yline(0) xline(0) border b2(" ");

set more 0; more ; set more 1 ;

di _n in wh _dup(75) "-" _n in gr
"Explanation: The rectangular impression is caused by the effects of a few" _n
"extreme outliers. Most likely, the outlier(s) in the x-direction will be " _n
"for different observations than the outliers in the y-direction. So " _n
"outliers cause big steps that are very nearly parallel either to the" _n 
"x-axis or to the y-axis."_n(2) 
"Here is another Cauchy sample ..." _n 
in wh _dup(75) "-" _n ;

di in wh ". replace sumC  = sum(invnorm(uniform())/invnorm(uniform()))" ;
push replace sumC  = sum(invnorm(uniform())/invnorm(uniform())) ;
nois replace sumC  = sum(invnorm(uniform())/invnorm(uniform())) ;

di in wh ". replace sumC2 = sum(invnorm(uniform())/invnorm(uniform()))" ;
push replace sumC2 = sum(invnorm(uniform())/invnorm(uniform())) ;
nois replace sumC2 = sum(invnorm(uniform())/invnorm(uniform())) ;

di in wh 
". graph sumC sumC2, c(l) s(i) ti(...) yline(0) xline(0) border b2( )" ;
push graph sumC sumC2, c(l) s(i) ti("Two dimensional Cauchy random walk") 
		       yline(0) xline(0) border b2(" ") ;
nois graph sumC sumC2, c(l) s(i) ti("Two dimensional Cauchy random walk") 
		       yline(0) xline(0) border b2(" ") ;

set more 0; more ; set more 1 ;

di _n in wh _dup(75) "-" _n in gr
"That is all. " _n 
in wh _dup(75) "-" _n ;

exit ;

First version of Tutorial was made by Albert Verbeek, June '90, based on 
Huber's ISP random walk demo. Modification by Jeroen Weesie (Apr '97)