How to Implement Newey-West covariance matrix properly for MDE estimation

I am trying to implement the MDE method for GARCH given by Baillie and Chung '01 (Estimation of GARCH Models from the Autocorrelations of the Squares of a Process, Jrl. of Time Series Analysis) but I can't understand it enough to implement the algorithm.

For a squared process $$X \sim GARCH(\theta)$$ with size $$n$$, MDE criterion is

$$\arg\min_\theta \left(R C^{-1}_{NW}R'\right)\tag1$$

where R equals to $$\hat\rho - \rho(\theta)$$ with $$\hat\rho$$ as sample autocorrelations and $$\rho(\theta)$$ as parametric autocorrelation function for GARCH process. And the $$C_{NW}$$ is the Newey-West Covariance matrix:

$$C_{NW} = \gamma_0^{-2} S\tag2$$

Here

$$S=\Omega_0 + \sum_{i=1}^m \left(1- \frac{i}{1+m}\right) \left(\Omega_i+\Omega_i'\right)\tag3$$

with $$\Omega_i= \frac{\sum_{t=i+1}^n Z_t Z_{t-i}'}{n}\tag4$$

And for $$D=x-\bar{x}$$

$$Z_t = \begin{bmatrix} D_t \, D_{t-1} - \rho_1(\theta)\,D_t^2\\ \vdots\\ D_t \, D_{t-k} - \rho_k(\theta)\,D_t^2 \end{bmatrix}\tag5$$

The value $$m$$ is the number of positive autocorrelations and for simplicity lets say $$k=5$$. $$Z_t$$ and $$\rho$$ are $$k\times 1$$ vectors, $$\Omega$$ and $$S$$ are $$k\times k$$ matrices.

I implemented the Eq.5 as follows

#--------------------------------------------------------
# for Eq. 5
# as a (k x k) Z matrix containing Z_1 to Z_k
z = np.zeros((k, k), np.float64)
t=0
while t<k:
i = t+1
z[t] = D[t]*D[i:k+i]-rho[:k]*D[t]**2
t = i
#----------------------------------------------------


I am not sure if this is correct way to build $$Z$$ (mathematically it is correct for $$Z_t$$, but what about $$Z_{t-i}$$?), and since $$k$$, $$m$$ and $$n$$ are not equal, I am confused about Eq. 3 and 4. So can someone help me and clearify things?

• Hi: I only slightly more than glanced at the paper but $k$ is the dimension of each $Z_t$ observation, $m$ is the number of lags you want to use when calculating the newey west estimator ( this is subjective. see the original NW paper for details ) and $n$ is the number of $Z_t$ observations in the sample. So, one would not necessarily expect the value of these parameters to be the same. Jul 16, 2023 at 7:18
• This paper below has some suggestions for how to deal with the choice of $m$. Notice that your problem is not totally standard Newey West because, usually, NW uses autocorrelations ($\rho_j$ are scalars in the paper below ) for the NW estimate. In your case, you have COVARIANCE MATRICES that need to be used for the estimate. So, Baillie and Chung are kind of doing an extension of Newey West or, if you don't want to call it an extension, it's an application of Newey West to the multivariate case. Jul 16, 2023 at 7:23
• HI Erdem: Three things. 1) I'll try to read the whole paper during this week because it really needs to be read completely and carefully. 2) I'm not sure how long 1) will take me ? 3) Even when I read it, I don't know if I'll be able to say anything helpful ? Maybe someone else on this list who works in this area can look at it and say something helpful ? I haven't worked with NW so thanks to anyone else who can help. I'll get back to you either way. Jul 16, 2023 at 23:08
• Hi Erdem: yes, the bartlett sample autocorrelation formula that looks confusing is actually not all that confusing when you apply it so specific models. It only looks confusing because of the summation from 1 to infinity. but, for specific models, the sum is a lot less. I'm still reading but I'm glad you are making some progress. I always wanted to learn more about GMM concept and using autocorrelations for robustness also so I'm going to keep going at my snail's pace. Jul 21, 2023 at 14:46
• @Bob Jansen let's see. I first tried to ask this question on this community but it was closed (the reason was "it is about programming"). As i was new in the community i thought it is not allowed, so i reedit a bit and searched another related community. I found finance... I waited an answer for 2 weeks and i researched at the same time to solve it. then i answered. later i learnt that it was my misunderstanding and it could be published here with some small edits. i have no knowledge on finance and i see my self more a statistician... so, here, CV is home to this question. Aug 26, 2023 at 9:14

I found the solution using a heuristic way. As mentioned in the comments obove, the size of the sum in ($$4$$) seems a bit problematic. After analysing the results, I took the upper bound of the sum as $$n-k$$ (Since $$Z_t$$ is $$k\times1$$, using ($$4$$) as it is, returns errors)

For the size of Bartlett's weights ($$m$$), i used $$m=k$$ (It was a personal choice, I decided by comparing this with the results of the formulas suggested in the literature)

The working codes of the implementation for GARCH process is shared below:

# MDE Newey-West
#----------------------------------------------------
#----------------------------------------------------
k=30                                        # Zt-size
k1= k+1
m = 30                   # size of Bartlett's weights
m1= m+1
#----------------------------------------------------
_x = x[::-1]               # for timeserries indexing
#----------------------------------------------------
scor = acf(_x, nlags=nobs)  # sample autocorrelations
#----------------------------------------------------
D = _x - x.mean()
#----------------------------------------------------
w = 1 - np.arange(m1)/(m1)       # Bartlett's weights
#====================================================
@jit
def MDE_NW(param) -> float:
#------------------------------------------------
a, b = param
atb = a+b
ab = a*b
#------------------------------------------------
if (a <= 0.0) or (b <= 0.0) or (1-atb <= 0):
return 1e12
#------------------------------------------------
# parametric autocorrelations
ro = np.zeros(k1)
ro[1] = a + (ab*a) / (1 - 2*ab - b**2)
i=2
while i<m:
ro[i] = ro[1] * (atb)**(i-1)
i+=1
pcor = ro[1:]
#------------------------------------------------
# Z matrix (here as h)
h = np.zeros((nobs-k, k))
t=0
while t<nobs-k:
i = t+1
h[t] = D[t]*D[i:k+i]-pcor*D[t]**2
t = i
#------------------------------------------------
R  = scor[1:k1]-pcor
#------------------------------------------------
# https://www.statsmodels.org/0.6.1/_modules/statsmodels/stats/sandwich_covariance.html#:~:text=S_hac_simple
S = dot(h.T, h)
i=1
while i<m1:   # m+1
omg = dot(h[i:].T, h[:-i])
S += w[i] * (omg + omg.T)
i+=1
return R.T @ inv(S) @ R


One can implement this function using scipy.optimize.minimize (BFGS and Nelder-Mead algorithms are suggested. BFGS is faster and mostly better but sometimes returns wrong values)

• Thanks Erdem. I'll check it out and hopefully learn something. I'm still going through the paper ( sllllllllllllooooooooooooooooooooooowly ). Jul 28, 2023 at 20:18