Purpose
The way to generate distribution for bias and Pseudo t statistic

Bias simulation

> set.seed(1977)
> results <- data.frame()
> i <- 1
> for (i in 1:5000) {
+     T <- 1000
+     inc <- rnorm(T)
+     y <- cumsum(inc)
+     num <- T^-1 * sum(y[1:(T - 1)] * inc[2:T])
+     denom <- T^-2 * sum(y[1:(T - 1)]^2)
+     temp <- num/denom
+     results <- rbind(results, temp)
+ }
> quantile(results[, 1], probs = c(0.01, 0.05, 0.1))
        1%         5%        10%
-14.371762  -8.183875  -5.866235

The above quantiles are same as the ones given in Fuller Table.

> par(mfrow = c(1, 2))
> hist(1000 * (results[, 1]), xlab = "", main = expression(paste(hat(delta),
+     "  Density")))
> plot(ecdf(1000 * (results[, 1])), xlab = "", main = expression(paste(hat(delta),
+     "  CDF")))

Bias-002.jpg

Pseudo t statistic Simulation

> set.seed(1977)
> results.t <- data.frame()
> i <- 1
> for (i in 1:5000) {
+     T <- 1000
+     inc <- rnorm(T)
+     y <- cumsum(inc)
+     num <- sum(y[1:(T - 1)] * inc[2:T])
+     denom <- sum(y[1:(T - 1)]^2)
+     temp <- num/denom
+     sigma.inc <- T^(-1) * (sum(inc^2))
+     sigma <- sigma.inc * sum(y[1:(T - 1)]^2)^(-1/2)
+     tau <- (temp)/sigma
+     results.t <- rbind(results.t, tau)
+ }
> quantile(results.t[, 1], probs = c(0.01, 0.05, 0.1))
       1%        5%       10%
-2.640077 -1.979768 -1.659552
> par(mfrow = c(1, 2))
> hist((results.t[, 1]), main = expression(paste(hat(tau), "  density")))
> plot(ecdf((results.t[, 1])), main = expression(paste(hat(tau),
+     "  CDF")))

Bias-004.jpg

Comparison of cdf Standard Normal Vs Tauu

> par(mfrow = c(1, 1))
> plot(ecdf(rnorm(T)), main = expression(paste(hat(tau), "  CDF")),
+     pch = 19, col.points = "red", xlim = c(-4, 4))
> par(new = T)
> plot(ecdf(results.t[, 1]), main = expression(paste(hat(tau),
+     "  CDF")), pch = 19, col.points = "blue", xlim = c(-4, 4))
> legend("topleft", legend = (c("StdNormal", "Tau")), fill = c("red",
+     "blue"))

Bias-005.jpg

Comparison of density Standard Normal Vs Tauu

> par(mfrow = c(1, 1))
> plot(density(rnorm(T)), main = "", col = "red", lwd = 4, xlim = c(-4,
+     4), ylim = c(0, 0.5), xlab = "", ylab = "")
> par(new = T)
> plot(density(results.t[, 1]), main = "", col = "blue", lwd = 4,
+     xlim = c(-4, 4), ylim = c(0, 0.5), xlab = "", ylab = "")
> legend("topleft", legend = (c("StdNormal", "Tau")), fill = c("red",
+     "blue"))

Bias-006.jpg