%% Estimate the Newey-West OLS Coefficient Covariance Matrix
% Model nominal GNP (GNPN) with consumer price index (CPI), real wages (WR), and the money stock (MS) using the linear model:
% \texttt{GNPN}_i = {\beta _0} + {\beta _1}\texttt{CPI}_i + {\beta _2}\texttt{WR}_i + {\beta _3}\texttt{MS}_i + \varepsilon_i.
% Estimate the model coefficients and the Newey-West OLS coefficient covariance matrix.
% Load the Nelson Plosser data set.

load Data_NelsonPlosser
Tbl = DataTable(:,[8,10,11,2]);  % Tabular array containing the variables
T = sum(~any(ismissing(Tbl),2)); % Remove NaNs to obtain sample size

y = Tbl{:,4};   % Numeric response
X = Tbl{:,1:3}; % Numeric matrix of predictors

%
% Fit the linear model. Remove the beginning block of NaN values in the residual vector for autocorr.

Mdl = fitlm(X,y);
resid = Mdl.Residuals.Raw(~isnan(Mdl.Residuals.Raw));
 
figure
subplot(2,1,1)
hold on
plotResiduals(Mdl,'fitted')
axis tight
plot([min(Mdl.Fitted) max(Mdl.Fitted)],[0 0],'k-')
title('Residual Plot')
xlabel('$\hat y$','Interpreter','latex')
ylabel('Residuals')
axis tight
subplot(2,1,2)
autocorr(resid)

% The residual plot exhibits signs of heteroscedasticity, autocorrelation, and possibly model misspecification. The sample autocorrelation function clearly exhibits autocorrelation.
%
% Calculate the lag selection parameter for the standard Newey-West HAC estimate (Andrews and Monohan, 1992).

maxLag = floor(4*(T/100)^(2/9));

%
% Estimate the standard Newey-West OLS coefficient covariance using hac by setting the bandwidth to maxLag + 1. Display the OLS coefficient estimates, their standard errors, and the covariance matrix.

EstCov = hac(X,y,'bandwidth',maxLag+1,'display','full');

% The first column in the output contains the OLS estimates ({\hat \beta _0},...,{\hat \beta _3}, respectively), and the second column contains their standard errors. The last four columns contained in the table represent the estimated coefficient covariance matrix. For example, Cov({\hat \beta _1},{\hat \beta _2}) =  - 0.2960.
%
% Alternatively, pass in a tabular array to hac.

EstCov = hac(Tbl,'bandwidth',maxLag+1,'display','off');

% The advantage of passing in a tabular array is that the top and left margins of the covariance table use the variable names.
% Copyright 2015 The MathWorks, Inc.