# Posted By

mjaniec on 10/26/11

# Statistics

Viewed 268 times
Favorited by 0 user(s)

# Drunk and her dog visualized

/ Published in: R
Save to your folder(s)

code inspired by the story presented in http://www-stat.wharton.upenn.edu/~steele/Courses/434/434Context/Co-integration/Murray93DrunkAndDog.pdf

Copy this code and paste it in your HTML
`library(tseries) 	# po.testlibrary(urca) 	# ca.jo Nmoves <- 1e5cfreq <- 0.01		      # correction frequencycfactor <- c(0.1,0.3,0.6)	# correction efficiency; 1,2 - length, 3 - angle; (0-1) ### drunk_path <- matrix(0,Nmoves,2)dog_path   <- matrix(0,Nmoves,2) random_walk <- rnorm(Nmoves*2,mean=0,sd=1) for (i in 2:Nmoves) {  	if (runif(1)>cfreq) { 		drunk_path[i,] <- drunk_path[i-1,]+rnorm(2,mean=0,sd=1)		dog_path[i,]   <- dog_path[i-1,]+rnorm(2,mean=0,sd=1) 	} 	else { 		d <- dog_path[i-1,]-drunk_path[i-1,] # delta 		h <- sqrt(d[1]^2+d[2]^2)		 # opposite	 		alpha <- atan2(d[2],d[1])		 # arc 		d1 <- h*runif(1,min=cfactor[1],max=cfactor[2])  	# reduced delta drunk		d2 <- h*runif(1,min=cfactor[1],max=cfactor[2])  	# reduced delta dog 		# a1 <- alpha*runif(1,min=cfactor[3],max=(2-cfactor[3]))	# distorted alpha drunk		# a2 <- alpha*runif(1,min=cfactor[3],max=(2-cfactor[3]))	# distorted alpha dog 		a1 <- alpha+(1-cfactor[3])*runif(1,-pi,pi)		a2 <- alpha+(1-cfactor[3])*runif(1,-pi,pi) 		d_drunk <- c(d1*cos(a1),d1*sin(a1))		d_dog   <- c(d2*cos(a2),d2*sin(a2)) 		drunk_path[i,] <- drunk_path[i-1,]+d_drunk		dog_path[i,]   <- dog_path[i-1,]-d_dog 		# cat("d=",d,"a=",alpha,a1,a2,"\n") 	} } par(mfrow=c(1,1)) xscope <- c(min(drunk_path[,1],dog_path[,1]),max(drunk_path[,1],dog_path[,1]))yscope <- c(min(drunk_path[,2],dog_path[,2]),max(drunk_path[,2],dog_path[,2])) plot(drunk_path[,1],drunk_path[,2],type="l",xlim=xscope,ylim=yscope,xlab="x",ylab="y")lines(dog_path[,1],dog_path[,2],col="Red") points(drunk_path[Nmoves,1],drunk_path[Nmoves,2],type="p")points(dog_path[Nmoves,1],dog_path[Nmoves,2],type="p") abline(v=drunk_path[Nmoves,1],col="Green")abline(h=drunk_path[Nmoves,2],col="Green")abline(v=dog_path[Nmoves,1],col="Blue")abline(h=dog_path[Nmoves,2],col="Blue") ### distance between pathes dd_path <- drunk_path-dog_path distance <- sqrt(dd_path[,1]^2+dd_path[,2]^2) plot(distance,type="l") ### mean(distance);sd(distance) drunk_path[Nmoves,]; dog_path[Nmoves,] delta <- drunk_path[Nmoves,]-dog_path[Nmoves,]sqrt(delta[1]^2+delta[2]^2) sum(drunk_path-dog_path)/Nmoves ### par(mfrow=c(4,1)) plot(drunk_path[,1],type="l") plot(drunk_path[,2],type="l") plot(dog_path[,1],type="l") plot(dog_path[,2],type="l") cor(drunk_path[,1],dog_path[,1]); cor(drunk_path[,2],dog_path[,2]); cor(drunk_path[,1],drunk_path[,2]) ### Phillips-Ouliaris Cointegration Test x_pathes <- cbind(drunk_path[,1],dog_path[,1])y_pathes <- cbind(drunk_path[,2],dog_path[,2])d_pathes <- cbind(drunk_path[,1],drunk_path[,2]) po.test(x_pathes)po.test(y_pathes)po.test(d_pathes) ### Johansen test colnames(x_pathes)<-c("drunk X","dog X")colnames(y_pathes)<-c("drunk Y","dog Y")colnames(d_pathes)<-c("drunk X","drunk Y") summary(ca.jo(x_pathes,type="eigen"))summary(ca.jo(y_pathes,type="eigen"))summary(ca.jo(d_pathes,type="eigen"))`