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$ 3 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} $$
that is a overdermined last squares problem, using $z = \Delta y$:
Ordinary Least Squares
$$ \mathbf{X} \cdot \mathbf{\beta} = z $$
The last squares solution is, :
$$ \mathbf{\hat{\beta}} = \left(\mathbf{X}^T \mathbf{X}\right)^{-1} \mathbf{X}^T z $$
The 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}} $$
The 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_t$ when I use Cholesky to solve $(X^T X)^{-1}$ I find some singular matrixes.
I was wondering if anything is wrong with my understanding, or some kind of regularization is recomended on $X$ or any kind of normalization in $y_t$ is requeried prior applying the test.
