# 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$$ 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.

Actually forming $$X^TX$$ and solving the linear system can lead to condition issues, since $$\kappa(A^TA)=\kappa(A)^2$$. Two numerically stable ways to solve the least squares problem $$A\beta = z$$ are by either using a QR decomposition or a singular value decomposition of $$A$$. Chances are that there are efficient implementations of these routines in your computin environment and you can see how to apply these matrix factorizations to solve the least squares problem here.

• thank you @whopowell96 I tried (Principal Component Regression)[ en.wikipedia.org/wiki/Principal_component_regression] on the least squares problem unfortunantly I found it prohibitive because SVD algorithm is not suitable for GPU optimization due it's nature, and made it hundred times slower. I am actually taking a look at QR factorization thank you I will see if it works.
– imbr
Feb 14, 2020 at 0:58

The errors vanished using a more numerically stable way like suggested by @whpowell96.

The following is the ADFs t statistics formulation using QR factorization

First perform factorization on $$X_{N_{eq},N_{p}}$$ matrix:

$$Q_{N_{eq},N_{p}} R_{N_{p},N_{p}} = X$$

For simplicity lets use $$N_{eq} = N$$ and $$N_{p} = P$$ here

$$Q_{N,P}R_{P,P} = X_{N, P}$$

Where $$R$$ is a square upper-triangular and $$Q$$ is orthogonal. The last squares solution is, :

$$\mathbf{\hat{\beta}} = \left(\mathbf{X}^T \mathbf{X}\right)^{-1} \mathbf{X}^T z$$ $$\mathbf{\hat{\beta}} = R^{-1} Q^T z$$

Sequence for numbericall stability:

1. Do the cross-product $$qtz = Q^T z$$
2. Use back-substitution to solve for the triangular $$R^{-1}$$ system for $$qtz$$
• Like with Pytorch torch.triangular_solve(qtz, R, upper=True)

The estimated $$\sigma^2$$ is the same for the classical OLS
The estimated standard error of each coeficient $$\mathbf{\hat{\beta}}_j$$ is

$$Var\left(\mathbf{\hat{\beta}}\right) = s^2 C$$

$$C = \left(X^T X\right)^{-1} = \left(R^T R\right)^{-1}$$

But for the standard error of the parameter $$j$$ we only care about jj'th diagonal entry of C. If $$e_j$$ is the j'th euclidean vector (all zeroes, except for a single one at the j'th position. You can write that entry as:

$$d_{jj} = e_j^T\cdot(X^TX)^{-1}\cdot e_j = e_j^T\cdot\left(R^T R\right)^{-1}\cdot e_j = (R^{-T}e_j)^T(R^{-T}e_j) = y^Ty$$

If you introduce the vector $$y = R^{-T}e_j =(R^{T})^{-1}e_j$$

To solve $$y$$ numerically use:

1. Form $$e_j$$
2. Use back-substitution to solve for the triangular $$R^{-T}$$ system for $$e_j$$
• Like with Pytorch torch.triangular_solve(ej, R.t(), upper=False)

Credits to here.

$$\hat{\sigma}_j = \hat{s.e.}\left(\mathbf{\hat{\beta}}_j \right) = \sqrt{s^2 C_{jj}} = \sqrt{s^2 d_{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 C_{33}} = \sqrt{s^2 d_{33}}$$

The the test statistic for $$\mathbf{\hat{\gamma}}$$ is

$$t = \frac{\mathbf{\hat{\gamma}}}{\hat{\sigma}_3}$$