Augmented Dickey Fuller (ADF) test statistics GPU formulation I have followed different sources of information and achieved the following formulation for the ADF $t$ test statistics. I implemented it to run several hundred thousands of ADFs $t$ statistics on GPU at once.
My formulation
Follows almost entirely Wikipedia convention:
$$ \Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \delta_1 \Delta y_{t-1} + ... + \delta_{p-1} \Delta y_{t-p+1} + \epsilon_t $$

*

*$ \epsilon_t $ error term (will be minimized)

*$\Delta$ first difference operator

*$\alpha$ drift term

*$\beta$ time trend or deterministic term

*$p$ lag order must be choosen, using for example the Akaike and Schwartz (Bayesian) information criteria.

For example for input data $y_t$ in {0, 1, 2 ... 8}, N=9 and p=2, the overdetermined system is:
$$ \Delta y_3 = \alpha + \beta 3 + \gamma y_{2} + \delta_1 \Delta y_2 + \delta_2 \Delta y_1 $$
$$ \Delta y_4 = \alpha + \beta 4 + \gamma y_{3} + \delta_1 \Delta y_3 + \delta_2 \Delta y_2  $$
$$ \Delta y_5 = \alpha + \beta 5 + \gamma y_{4} + \delta_1 \Delta y_4 + \delta_2 \Delta y_3  $$
$$ \Delta y_6 = \alpha + \beta 6 + \gamma y_{5} + \delta_1 \Delta y_5 + \delta_2 \Delta y_4  $$
$$ \Delta y_7 = \alpha + \beta 7 + \gamma y_{6} + \delta_1 \Delta y_6 + \delta_2 \Delta y_5  $$
$$ \Delta y_8 = \alpha + \beta 8 + \gamma y_{7} + \delta_1 \Delta y_7 + \delta_2 \Delta y_6  $$

*

*General case with data points $y_t$ in {0, 1, 2 ... N} and p=p.


*Parameters are $\alpha$, $\beta$, $\gamma$ plus p $\delta$'s, total is $3+p$.


*N must be $> 3+p$, otherwise cannot apply ordinary least squares.
$$
 \begin{pmatrix}
  1 &  p+1  & y_{p}    &  \Delta y_{p}    & \Delta y_{p-1}&  \cdots  &   \Delta y_{1} \\
  1 &  p+2  & y_{p+1}  &  \Delta y_{p+1}  & \Delta y_{p}  &  \cdots  &   \Delta y_{2} \\
  1 &  p+3  & y_{p+2}  &  \Delta y_{p+2}  & \Delta y_{p+1}&  \cdots  &   \Delta y_{3} \\ 
  \vdots &  \vdots  & \vdots  &  \vdots  &  \vdots  & \ddots  & \vdots  \\
  1 &  N-1  & y_{N-2}  &  \Delta y_{N-2}  & \Delta y_{N-3}&  \cdots  &   \Delta y_{N-p-1} \\ 
 \end{pmatrix} 
 \begin{pmatrix}
    \alpha   \\
    \beta     \\
    \gamma    \\
    \delta_1  \\
    \delta_2  \\
    \vdots \\
    \delta_p \\
\end{pmatrix} 
 = 
  \begin{pmatrix}
  \Delta y_{p+1}\\
  \Delta y_{p+2}\\
  \Delta y_{p+3}\\
  \vdots \\
  \Delta y_{N-1}\\
 \end{pmatrix}
$$
$$ \mathbf{X}_{N-p-1,3+p} \cdot \mathbf{\beta}_{3+p} =  \Delta y_{N-p-1}, $$
which is a overdetermined least squares problem, using $z = \Delta y$:
Ordinary Least Squares
$$ \mathbf{X} \cdot \mathbf{\beta} =  z $$
The least squares solution is:
$$ \mathbf{\hat{\beta}} = \left(\mathbf{X}^T \mathbf{X}\right)^{-1} \mathbf{X}^T  z $$
with error
$$   \mathbf{\hat{\epsilon}} =  z - \mathbf{X} \mathbf{\hat{\beta}}.$$
The number of linear equations is
$$ N_{eq} = N-p-1, $$
the number of parameters being estimated is ($\alpha$, $\beta$, $\gamma$ plus p $\delta$'s):
$$ N_p = (p+3). $$
The reduced chi-squared statistic, estimated $\sigma^2$ :
$$   \hat{\sigma}^2 =  s^2 = \frac{\mathbf{\hat{\epsilon}}^T \mathbf{\hat{\epsilon}}}{N_{eq}-N_p}, $$
the denominator is the statistical degrees of freedom.
If the strict exogeneity does not hold (as is the case with many time series models, where exogeneity is assumed only with respect to the past shocks but not the future ones), then these estimators will be biased in finite samples.
The estimated standard error of each coeficient  $\mathbf{\hat{\beta}}_j$ is
$$ Var\left(\mathbf{\hat{\beta}}\right)   = \sigma^2 Q = s^2 \left(\mathbf{X}^T \mathbf{X}\right)^{-1}.$$
The estimate of this standard error is obtained by replacing the unknown quantity $\sigma^2$ with its estimate $s^2$.
Thus the standard error for the parameter $j$ is:
$$ \hat{\sigma}_j = \hat{s.e.}\left(\mathbf{\hat{\beta}}_j \right)  = \sqrt{s^2 \left(\mathbf{X}^T \mathbf{X}\right)^{-1}_{jj}}.$$
For the specific case of $\mathbf{\hat{\beta}}_3$ ~ $\gamma$:
$$ \hat{\sigma}_3 = \hat{s.e.}\left(\mathbf{\hat{\gamma}} \right)  = \sqrt{s^2 \left(\mathbf{X}^T \mathbf{X}\right)^{-1}_{33}}.$$
Then the test statistic for $\mathbf{\hat{\gamma}}$ is:
$$ t = \frac{\mathbf{\hat{\gamma}}}{\hat{\sigma}_3}. $$
The problem
Unfortunately, for some specific data sets $y$ when I use Cholesky to solve $(X^T X)^{-1}$ I find some singular matrices.
I was wondering if anything is wrong with my understanding, or some kind of regularization is recommended on $X$ or any kind of normalization in $y$ is required prior to applying the test.
A piece of code Pytorch
X = th.zeros((batch, n-p-1, 3+p), device=dev, dtype=th.float32)
y = th.tensor(data, device=dev, dtype=th.float32)
diffilter = th.tensor([-1., 1.], device=dev, dtype=th.float32).reshape(1, 1, 2)
y = y.reshape(batch, 1, -1)
dy = th.conv1d(y, diffilter).reshape(batch, -1)
y = y.reshape(batch, -1)
z = dy[:, p:].clone()
if verbose:
    print(len(z), nobs, p, n, X.shape)
# X matrix
X[:, :, 0] = 1 # drift
X[:, :, 1] = th.arange(p+1, n) # deterministic trend
X[:, :, 2] = y[:, p:-1]# regression data (nobs)
# fill in columns, max lagged serial correlations
for i in range(1, p+1):
    X[:, :, 2+i] = dy[:, p-i:-i]

Xt = X.transpose(dim0=1, dim1=-1)
z = z.unsqueeze(2)
L = th.cholesky(th.bmm(Xt, X))
Gi =  th.cholesky_solve(th.eye(p+3), L) # ( X^T . X ) ^-1
Xtz = th.bmm(Xt, z)
Bhat = th.bmm(Gi, Xtz)
er = z - th.bmm(X, Bhat)
s2 = (th.matmul(er.transpose(dim0=1, dim1=-1), er)/(nobs-(p+3))).view(-1)
Bhat = Bhat.squeeze(2)
tstat = Bhat[:, 2]/th.sqrt(s2*Gi[:,2,2])

 A: One candidate reason is multicollinearity, which arises when some linear combination of your regressors yields another regressors. For example, it might happen that (although the chance is quite small, but then again if you run many replications, it is enough if it happens once)
$$
\begin{pmatrix}
\Delta y_{p}    & \Delta y_{p-1}&  \\
\Delta y_{p+1}  & \Delta y_{p}  \\
\Delta y_{p+2}  & \Delta y_{p+1}\\ 
 \vdots  & \vdots  \\
  \Delta y_{N-2}&  \Delta y_{N-3}
\end{pmatrix}
$$
has identical columns when $N$ is small, because the changes of the series happen to be constant over the rows.
A: Just to record the solution. @Christoph Hanck was right. The problem was with Pytorch/Libtorch not letting me discard failed cholesky. So I had to use the following simple cholesky for Pytorch. 
# Batched Cholesky nograd_cholesky
def cholesky(A):
    L = th.zeros_like(A)

    for i in range(A.shape[-1]):
        for j in range(i+1):
            s = 0.0
            for k in range(j):
                s = s + L[...,i,k] * L[...,j,k]

            L[...,i,j] = th.sqrt(A[...,i,i] - s) if (i == j) else \
                      (1.0 / L[...,j,j] * (A[...,i,j] - s))
    return L

