N <- 100000
p <- 0.5

Xj <- 0
X <- numeric(N)
for (j in 1:N) {
  V <- runif(1)
  Y <- ifelse(V <= p, Xj + 1, Xj - 1)
  U <- runif(1)

  if (Y > Xj) {
    alpha <- min(1, 2^(abs(Xj) - abs(Y))*(1-p)/p)
  } else {
    alpha <- min(1, 2^(abs(Xj) - abs(Y))*p/(1-p))
  }

  if (U <= alpha) {
    Xj <- Y
  }
  X[j] <- Xj
}

hist(X, breaks=seq(-20.5, 20.5, by=1), prob = TRUE)
xx <- seq(-15, 15, by=1)
points(xx, 1/3*2^(-abs(xx)), col="red")

plot(X[1:1000], type="l")