Simulate the critical values of Dickey-Fuller test statistic with R

Question

I am trying to replicate (approximately) the results of the unit root tests described in chapter 17 of Hamilton's book. For the case

\begin{equation} y_{t}=y_{t-1}+\varepsilon_{t},\ \varepsilon_{t}\sim \mathcal{N}(0,1) \end{equation}

with the statistic

\begin{equation} t=\frac{\hat{\rho}-1}{\hat{\sigma}_{\hat{\rho}}} \end{equation}

where $\hat{\rho}$ is the usual OLS estimate and $\hat{\sigma}_{\hat{\rho}}$ its standard error.

My question is:

• Is my R code correct? My results are not close to the tabulated values, for example

The critical values of the test with $T=25$ and $\alpha=0.05$ are -2.26 and 2.16. But I get -0.1433872 and 3.2984347

The function simulates $R$ times a random walk of size $T$ and calculates the least squares estimator for each, and their standard error. The resulting distribution should be similar to the tabulated values commonly used.

    case1<-function(nob,R){
      R<-1000
      beta.hat<-statistic<-sd.beta<-c()
      for (r in 1:R){
          y<-error<-vector()
          y[1]<-0
          nob<-25
          for (i in 2:nob){
            y[i]<-y[i-1]+rnorm(1,0,1)
          }
          res<-lm(y[2:nob]~Hmisc::Lag(y)[2:nob])
          res.sum <- summary(res)
          statistic[r] <- (res.sum$coefficients[2, 1]-1)/res.sum$coefficients[2,2]
      }
      statistic
    }
    data<-case1(nob=25,R=10000)
    quantile(data, probs=c(0.025,0.975))
    hist(data)

Thanks in advance.