## Posted By

mjaniec on 10/26/11

# Drunk and her dog visualized

/ Published in: R

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

1. library(tseries) # po.test
2. library(urca) # ca.jo
3.
4. Nmoves <- 1e5
5. cfreq <- 0.01 # correction frequency
6. cfactor <- c(0.1,0.3,0.6) # correction efficiency; 1,2 - length, 3 - angle; (0-1)
7.
8. ###
9.
10. drunk_path <- matrix(0,Nmoves,2)
11. dog_path <- matrix(0,Nmoves,2)
12.
13. random_walk <- rnorm(Nmoves*2,mean=0,sd=1)
14.
15. for (i in 2:Nmoves) {
16.
17. if (runif(1)>cfreq) {
18.
19. drunk_path[i,] <- drunk_path[i-1,]+rnorm(2,mean=0,sd=1)
20. dog_path[i,] <- dog_path[i-1,]+rnorm(2,mean=0,sd=1)
21.
22. }
23.
24. else {
25.
26. d <- dog_path[i-1,]-drunk_path[i-1,] # delta
27.
28. h <- sqrt(d[1]^2+d[2]^2) # opposite
29.
30. alpha <- atan2(d[2],d[1]) # arc
31.
32. d1 <- h*runif(1,min=cfactor[1],max=cfactor[2]) # reduced delta drunk
33. d2 <- h*runif(1,min=cfactor[1],max=cfactor[2]) # reduced delta dog
34.
35. # a1 <- alpha*runif(1,min=cfactor[3],max=(2-cfactor[3])) # distorted alpha drunk
36. # a2 <- alpha*runif(1,min=cfactor[3],max=(2-cfactor[3])) # distorted alpha dog
37.
38. a1 <- alpha+(1-cfactor[3])*runif(1,-pi,pi)
39. a2 <- alpha+(1-cfactor[3])*runif(1,-pi,pi)
40.
41. d_drunk <- c(d1*cos(a1),d1*sin(a1))
42. d_dog <- c(d2*cos(a2),d2*sin(a2))
43.
44. drunk_path[i,] <- drunk_path[i-1,]+d_drunk
45. dog_path[i,] <- dog_path[i-1,]-d_dog
46.
47. # cat("d=",d,"a=",alpha,a1,a2,"\n")
48.
49. }
50.
51. }
52.
53. par(mfrow=c(1,1))
54.
55. xscope <- c(min(drunk_path[,1],dog_path[,1]),max(drunk_path[,1],dog_path[,1]))
56. yscope <- c(min(drunk_path[,2],dog_path[,2]),max(drunk_path[,2],dog_path[,2]))
57.
58. plot(drunk_path[,1],drunk_path[,2],type="l",xlim=xscope,ylim=yscope,xlab="x",ylab="y")
59. lines(dog_path[,1],dog_path[,2],col="Red")
60.
61. points(drunk_path[Nmoves,1],drunk_path[Nmoves,2],type="p")
62. points(dog_path[Nmoves,1],dog_path[Nmoves,2],type="p")
63.
64. abline(v=drunk_path[Nmoves,1],col="Green")
65. abline(h=drunk_path[Nmoves,2],col="Green")
66. abline(v=dog_path[Nmoves,1],col="Blue")
67. abline(h=dog_path[Nmoves,2],col="Blue")
68.
69. ### distance between pathes
70.
71. dd_path <- drunk_path-dog_path
72.
73. distance <- sqrt(dd_path[,1]^2+dd_path[,2]^2)
74.
75. plot(distance,type="l")
76.
77. ###
78.
79. mean(distance);sd(distance)
80.
81. drunk_path[Nmoves,]; dog_path[Nmoves,]
82.
83. delta <- drunk_path[Nmoves,]-dog_path[Nmoves,]
84. sqrt(delta[1]^2+delta[2]^2)
85.
86. sum(drunk_path-dog_path)/Nmoves
87.
88. ###
89.
90. par(mfrow=c(4,1))
91.
92. plot(drunk_path[,1],type="l")
93.
94. plot(drunk_path[,2],type="l")
95.
96. plot(dog_path[,1],type="l")
97.
98. plot(dog_path[,2],type="l")
99.
100. cor(drunk_path[,1],dog_path[,1]); cor(drunk_path[,2],dog_path[,2]); cor(drunk_path[,1],drunk_path[,2])
101.
102. ### Phillips-Ouliaris Cointegration Test
103.
104. x_pathes <- cbind(drunk_path[,1],dog_path[,1])
105. y_pathes <- cbind(drunk_path[,2],dog_path[,2])
106. d_pathes <- cbind(drunk_path[,1],drunk_path[,2])
107.
108. po.test(x_pathes)
109. po.test(y_pathes)
110. po.test(d_pathes)
111.
112. ### Johansen test
113.
114. colnames(x_pathes)<-c("drunk X","dog X")
115. colnames(y_pathes)<-c("drunk Y","dog Y")
116. colnames(d_pathes)<-c("drunk X","drunk Y")
117.
118. summary(ca.jo(x_pathes,type="eigen"))
119. summary(ca.jo(y_pathes,type="eigen"))
120. summary(ca.jo(d_pathes,type="eigen"))