library(tseries) 	# po.test
library(urca) 	# ca.jo

Nmoves <- 1e5
cfreq <- 0.01		      # correction frequency
cfactor <- 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"))

Nmoves <- 1e5
cfreq <- 0.05		      # correction frequency
cfactor <- c(0.2,0.5,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")

	}

}

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(h=0,col="Green")
# abline(v=0,col="Green")
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")


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

http://www.reakkt.com/2011/10/cointegrated-drunk-and-her-dog.html

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

Drunk and her dog visualized

R