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

Question

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?

Answer 1

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)