Assume that I have a process $y_t$ such that $$y_t = c + \phi_1 y_{t-1} + \ldots + \phi_p y_{t-p} + u_t$$ where $u_t$ is i.i.d. white noise such that $E[u_t] = 0, \forall t$ and $E[u_t u_s]$ is equal to $\sigma^2$ if $t = s$ and 0 otherwise. Define $$x_t = \begin{bmatrix} 1 \\ y_{t-1} \\ \vdots \\ y_{t-p} \end{bmatrix} \quad \text{and} \quad \phi = \begin{bmatrix} c \\ \phi_{1} \\ \vdots \\ \phi_{p} \end{bmatrix}$$ so that we can write $$y_t = x_t'\phi + u_t.$$ This process can be estimated by OLS such that $$\hat{\phi} = \left(\sum_{t = 1}^{T} x_t x_t' \right)^{-1} \sum_{t = 1}^{T} x_t y_t$$ and $\sigma^2$ can be estimated by $$s^2 = \frac{1}{T - p - 1} \sum_{t = 1}^{T} \left(y - x_t'\hat{\phi} \right)^2.$$
Note that this analysis requires that we have access to $T + p$ samples of $y_t$ denoted as $(y_{-p+1}, y_{-p+2}, \ldots, y_0, y_1, \ldots, y_T)$. However, I have noticed that implementation of this result is very ambiguous. In a practical setting, you would have only $T$ observations from which estimation would be carried out using $T-p$ so that all the previous summations go from $p+1$ to $T$ and thus the correction should be done using $T - 2p - 1$. Using the R command ar with method ols uses $T - p$, while in a multivariate setting, package vars uses $T - kp - 1$ where $k$ is the number of endogenous variables.
Could someone explain the intuition behind these corrections over what I would assume to be the correct one (using $T - 2p - 1$)? I know that asymptotically it is all the same, but the corrections are done precisely in order to improve finite-sample properties and so I'm not sure why you wouldn't use such results.
