Estimating Covariance Matrix of Innovations of Multivariate Random Walk Suppose that I have a multivariate random walk:
$X_{t+1} = X_t + \epsilon_t$ where $\epsilon_t \sim N(0,\Sigma)$
Estimating the covariance matrix $\Sigma$ is straightforward from first differences $X_{t+1}-X_t$ using the sample covariance.
Instead if I would like to estimate the covariance directly from the levels $X_0,...,X_n$ again using the sample covariance but this time on the levels and then using the relation $Cov(X_n-X_0) = n \Sigma$ to find $\Sigma$ by simply dividing by $n$, I have a feeling that I will underestimate the variances hence overestimate the correlations. I also cannot decide whether to subtract the unconditional mean $X_0$ or sample mean of the levels during the calculation of covariance matrix. 
Basically does $\hat{Cov}(X_n-X_0) = \frac{1}{N} \Sigma_1^N(X_n-X_0)(X_n-X_0)^T$ or the same formula where $X_0$ replaced with sample average $\hat{\mu} = \frac{1}{N} \Sigma_1^N X_i$ makes more sense (in this scenario we divide by $N-1$ because of degrees of freedom but that's not the issue)?
Would I be underestimating the variances and overestimating the correlations?
Obviously my intention is not to estimate the covariance of innovations in this manner where there is already a way to estimate them using the difference series. My main question is how this is affecting OLS estimates of multiple regression when regressors have unit roots, as 
$\beta = \large{\frac{Cov(X,y)}{Cov(X)}}$
doesn't seem to be so well defined anymore? Yes they are super consistent in case of a cointegration, but that's a little after the fact in terms of estimation. 
 A: Consider the following single-equation model:
$$y_t = \rho y_{t-1} + \beta x_t + u_t,\;\;\; t=0,...,T,\;\;\;y_0=0 \tag{1}$$
with the following stochastic assumptions:
1) The $x$'s come from an i.i.d. sample
2) The $x$'s are strictly exogenous to the disturbance
3) The disturbance $u_t$ is white noise.
4) It follows but let's add it for clarity, that current $x$'s (sometimes called "forcing factors" in the ARMA-X literature) are independent of past dependent variables also.
Denote also $M_x = I - X(X'X)^{-1}X'$, with appropriate dimensions to fit the lagged structure, which is the residual maker matrix associated with the regressor $X$, a symmetric and idempotent matrix.
On the side, denote $S_t = \sum_{i=1}^t x_i,\;\;\; U_t = \sum_{i=1}^t u_i$. Then, if in reality $\rho =1$ (a unit root), we could re-write the model for $t-1$ as
$$y_{t-1} = \beta S_{t-1} + U_{t-1} \tag{2}$$
Keep that.
By partitioned-regression results (which are algebraic and do not depend on stochastic assumptions), we have that
$$\hat \rho_{OLS} = \left(\mathbf y_{-1}'M_x\mathbf y_{-1}\right)^{-1}\mathbf y_{-1}'M_x\mathbf y = \left((M_x\mathbf y_{-1})'(M_x\mathbf y_{-1})\right)^{-1}(M_x\mathbf y_{-1})'M_x\mathbf y \tag{3}$$
..the second equality due to the symmetry and idempotency of $M_x$.
Note that $\mathbf y_{-1}'M_x$ are the residuals (in row-vector form) from regressing $y_{t-1}$ on ${x_t}$, in an OLS regression $y_{t-1}= ax_t + v_{t-1}$. But given our stochastic assumptions, and as it should be obvious from $(2)%$, we have $a=0$, and the OLS regression will have no problem in detecting that. Which means that the residuals from this regression will tend to equal the dependent variable : $\mathbf y_{-1}'M_x \to \mathbf y_{-1}$.  
It follows that
$$\hat \rho_{OLS} \to \left(\mathbf y_{-1}'\mathbf y_{-1}\right)^{-1}\mathbf y_{-1}'(M_x\mathbf y) \tag{4}$$
Moreover, $M_x\mathbf y$ are the residuals from regressing $y_t$ on $x_t$, and, given our stochastic assumptions it is easy to conclude that $M_x\mathbf y \to \rho \mathbf y_{-1} + \mathbf u$. 
With this we obtain
$$\hat \rho \to \rho +  \left(\mathbf y_{-1}'\mathbf y_{-1}\right)^{-1}\mathbf y_{-1}'\mathbf u \to \rho \tag{5}$$
...because given our stochastic assumptions, current distrubances are orthogonal to past dependent variables. So we see that, irrespective of whether $\rho=1$ or not, the OLS estimator will estimate it consistently.
What about the $\beta$ coefficient? Aplying again partitioned regression results, and analogous notation as before, we have that 
$$\hat \beta_{OLS} = \left(X'M_{-1}X\right)^{-1}X'M_{-1}\mathbf y \tag{6}$$
Using previous reasoning $X'M_{-1} \to X'$. Using the specification, we have 
$$\hat \beta_{OLS} \to \left(X'X\right)^{-1}X'\mathbf y = \rho\left(X'X\right)^{-1}X'\mathbf y_{-1} + \beta + \left(X'X\right)^{-1}X'\mathbf u$$
With the same steps and reasoning as before we get that the first and third term will go to zero (always given our stochastic assumptions), and so 
$$\hat \beta_{OLS} \to \beta \tag{7}$$
I have not touched on the matter regarding the distribution of the estimator.
I will return as regards the rate of convergence and wether the estimator is superconsistent or not (and if it is for which coefficients).
