Simulation Dickey-Fuller distribution I tried to simulate the Dickey- Fuller statistics with a short Matlab routine. 
I set Y0=0 and generated
$$
Y_{t}=Y_{t-1} + \varepsilon_t,
$$
where $\epsilon_t$ is  $\mathcal{N}(0,1)$-distributed. Then I estimate the model
$$
\Delta Y_t = \beta Y_{t-1} + \varepsilon_t
$$
and computed the Dickey-Fuller statistic in order to estimate the quantile.



 function [Y] = genY(n)
 %generate random walk.
 Y=zeros(1,n);
 Y(1)=0;
 for i=2:n
     Y(i)=Y(i-1)+randn;
 end
 end




 function [dist] = estBeta(n)
 %estimate  Dickey-Fuller statistics in one step
 Y=genY(n);
 DeltaY=Y(2:n)-Y(1:(n-1));
 est=(sum(DeltaY(2:end).*DeltaY(1:(end-1))))/(sum((DeltaY(1:(end-1))).^2));
 dist=est/(std(DeltaY)/sqrt(n-1));
 end 




 function [quantil] = perc(m)
 %calculation of quantile
 z=zeros(1,m);
 for i=1:m
    z(i)=estBeta(1000);
 end
 quantil=quantile(z,0.01);
 end


If I execute the program I will always end up with -2.3 at p=0.01.
Has anybody an idea how to fix this problem?
I have a solution for the no contant no trend case (you can adapt the code for further model specifications). 


Set up the number of replications, I have done this in a for loop so that you can change it to meet your requirements. 
Deriving the DF tests state see W.

%% Replicating Dickey and Fuller critical values


for q = [25,50,100,250,500]
    T = q; 
    reps = 5000000;
    DFstats = zeros(1,reps);

for i = 1:reps
    u = normrnd(0,1,[1,T]);
    W = 1/sqrt(T)*cumsum(u);
    DFstats(i)=W(T)^2-1/(2*sqrt(mean(W.^2)));
end 

CV = sort(DFstats);
ActualValues = CV([(0.01*reps),(0.05*reps),(0.1*reps)]);
ActualTable(:,q)=ActualValues;


end
I have a solution for the no contant no trend case (you can adapt the code for further model specifications). 


Set up the number of replications, I have done this in a for loop so that you can change it to meet your requirements. 
Deriving the DF tests state see W.

%% Replicating Dickey and Fuller critical values


for q = [25,50,100,250,500]
    T = q; 
    reps = 5000000;
    DFstats = zeros(1,reps);

for i = 1:reps
    u = normrnd(0,1,[1,T]);
    W = 1/sqrt(T)*cumsum(u);
    DFstats(i)=W(T)^2-1/(2*sqrt(mean(W.^2)));
end 

CV = sort(DFstats);
ActualValues = CV([(0.01*reps),(0.05*reps),(0.1*reps)]);
ActualTable(:,q)=ActualValues;


end