Estimate HAC Covariance Matrix from data by hand - Newey West Given $T$ realizations for $N$ random variables, $X\in\mathbb{R}^{N\times T}$, I want to estimate the covariance matrix of the data, $\Omega\in\mathbb{R}^{N\times N}$. The sample covariance would be $$\hat\Omega_{i,j}=\frac{1}{T-1}\sum_{t=1}^T (X_{i,t}-\bar{X_i})(X_{j,t}-\bar{X_j}).$$
I however want to compute the hac (heteroscedasticity and autocorrelation) robust covariance matrix, using for example a Bartlett kernel. This is suggested by Newey and West (1987). Their paper is quite technical though and I am not sure how to compute the components of the hac covariance matrix. What would be the formula for $\hat\Omega_{i,j}$ in this case?

According to this description,
\begin{align}
\hat\Omega = \frac{T}{T-K}\sum_{j=-\infty}^\infty k(j/b_T)\hat\Gamma_j, \tag{1}
\end{align}
where

*

*$k$ is the kernel function with bandwitdth $b_T>0$,

*$K$ is the number of parameters and

*$\hat\Gamma_j$ are sample autocorrelations,
\begin{align*}
\hat\Gamma_j=\begin{cases}\frac{1}{T}\sum\limits_{t=1+j}^TX_tX_{t-j}' && j\geq0 \\ \\
\frac{1}{T}\sum\limits_{t=1-j}^TX_t'X_{t+j} && j<0. \end{cases}
\end{align*}
These are equations (F.4) and (F.5) from here. Are these formuale correct? Can I set $K=1$ like in the sample covariance matrix? For what $j$ do I truncate the sum in (1)?
 A: This is probably coming in late, but here's the basic idea.
You have some possibly vector-valued time series $\{Z_t\}_{t=1}^T$. More often than not, those would be score functions. For example, if you worked with a linear regression model $Y = X\beta + \epsilon$ (observations as rows), then  you'd use $Z_t = X_t \epsilon_t$ where replacing error terms with estimated residuals makes the estimators we'll see below feasible. Whether this is the case or not, you're interested in estimating the long-run covariance matrix of your vector-valued time series $\{Z_t\}$, i.e.
$$ \Omega_T := \underset{T \rightarrow \infty} \lim T \; \mathbb{E}\left( \bar{Z} \bar{Z}'  \right) $$
where I use $\bar Z := \sum_{t=1}^T Z_t / T$. (The $T$ appears outside because you're normalizing by $\sqrt{T}$). So, how do you do that? Typically, you'll use a weighted sum of sample autocovariances
$$ \Phi_T(\tau) := \sum_{t=\tau+1}^T Z_t Z_{t-\tau}' \; \; 0 \leq \tau \leq T-1 $$
Obviously, you're estimating a covariance matrix so the time series is covariance stationary and we have $\Phi_T(\tau) = \Phi_T(-\tau)'$ for $\tau < 0$. Now, you've probably realized $\Phi_T(T-1)$ would be estimated with just one observation which would most likely not be very precise, so what we tend to do is select weights based on so-called kernels which decay as $\tau$ increases. For Newey and West (1987), for a truncation lag of order $m_T > 0$, the weights are
$$ w(\tau) = 1 - \frac{|\tau|}{m_T + 1}. $$
So, we just need to apply our definition of long-run covariance and replace everything with weighted sample counterparts:
\begin{align}
    \hat{\Omega}_T &= \sum_{\tau=-m_T}^{m_T} \left( 1 - \frac{|\tau|}{m_T + 1} \right) \Phi_T(\tau) \\
                   &= \Phi_T(0) + \sum_{\tau=1}^{m_T} \left( 1 - \frac{|\tau|}{m_T + 1} \right) \left( \Phi_T(\tau) + \Phi_T(\tau)' \right)
\end{align}
where I exploited the symmetries $\Phi_T(\tau) = \Phi_T(-\tau)'$ and $w(\tau) = w(-\tau)$. This gives you directly a simple way to code this. For example, if you have observations as rows and variables as column in an object 'Z' in MATLAB, then
Omega = Z'.Z/Tz;
for jj = 1:mT

    % Compute sample autocovariance of order jj
    % NB: Need Zt and Zt-jj to be of the same dim.
    Z0   = Z((1+jj):end,:); 
    Zlag = Z(1:(end-jj),:);
    PhiT = Z0'*Zlag/size(Z0,1);
    
    % Weights (Bartlett)
    wj = 1 - jj/(mT + 1);

    % Update long-run covariance matrix estimate
    Omega = Omega + wj*( PhiT + PhiT' );
end

Basically, you start with $\Phi_T(0)$ and then you add up the weighted autocovariances of order 1 through $m_T$. To choose the truncation lag, you can use a rule of thumb like the nearest integer to $0.75T^{1/3}$. You can do more complicated things (optimal bandwidth selection would then require proceeding in steps because the formulas usually involve the unknown parameters you're trying to estimate here).
