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])
  • $\begingroup$ Can you report a seed and a code snippet to reproduce the problem? Also, can you elaborate on why the GPU part of the question is relevant? Would the problem not also occur, at some point, on a standard CPU? $\endgroup$ Feb 13, 2020 at 15:27
  • $\begingroup$ @ChristophHanck thank you. The GPU part is relevant like I mention bellow due not all possible algorithms for solving the least squares problem being GPU efficient. I will get an piece of code to paste here. $\endgroup$
    – iambr
    Feb 13, 2020 at 17:30

2 Answers 2


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.

  • $\begingroup$ Thanks @Christoph Hanck I thought that theoretically was possible too. Any ideas on how to solve that? I tried (Principal Component Regression)[ en.wikipedia.org/wiki/Principal_component_regression] on the least squares problem unfortunantly I found it prohibitive because SVD is not suitable for GPU optimization due it's nature, and made it hundred times slower. $\endgroup$
    – iambr
    Feb 13, 2020 at 17:26
  • 1
    $\begingroup$ Since your implementation is in principle correct and the issue is just "bad luck" drawing an unusual sample, I would just catch and then discard these cases. $\endgroup$ Feb 13, 2020 at 18:32
  • $\begingroup$ Unfortunantly I cannot discard any of the cases. And if I use QR factorization just figure out that problem doesn't happen. I was thinking on contaminate my data y_t with 1% of random noise to see if works. $\endgroup$
    – iambr
    Feb 14, 2020 at 13:47
  • $\begingroup$ Can you elaborate why you cannot discard cases? Jittering your data may be OK, but may also affect the statistical properties of the procedure, so that would not be my first choice. $\endgroup$ Feb 18, 2020 at 9:37
  • $\begingroup$ Yes you are right I can discard but I am not being able for now using pytorch GPU I will try that path a little more before giving up. Thanks again! $\endgroup$
    – iambr
    Feb 18, 2020 at 13:31

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.