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
  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"))

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

Report this snippet


Comments

RSS Icon Subscribe to comments

You need to login to post a comment.