%% Assess Conditional Heteroscedasticity Using the Lagrange Multiplier Test
% Test whether there are significant ARCH effects in a simulated response series using lmtest. The parameter values in this example are arbitrary.
% Specify the AR(1) model with an ARCH(1) variance:
% y_t = 0.9y_{t-1}+\varepsilon_t,
% where
% \varepsilon_t = w_t\sqrt{h_t}.
% h_t = 1 + 0.5\varepsilon^2_{t-1}.
% w_t is Gaussian with mean 0 and variance 1.

VarMdl = garch('ARCH',0.5,'Constant',1);
Mdl = arima('Constant',0,'Variance',VarMdl,'AR',0.9);

% Mdl is a fully specified, AR(1) model with an ARCH(1) variance.
%
% Simulate presample and effective sample responses from Mdl.

T = 100;
rng(1);  % For reproducibility
n = 2;   % Number of presample observations required for the gradient
[y,ep,v] = simulate(Mdl,T + n);

% ep is the random path of innovations from VarMdl. The software filters ep through Mdl to yield the random response path y.
%
% Specify the restricted model and assume that the AR model constant is 0:
% y_t = c + \phi_1y_{t-1}+\varepsilon_t,
% where h_t = \alpha_0 + \alpha_1\varepsilon^2_{t-1}.

VarMdl0 = garch(0,1);
VarMdl0.ARCH{1} = 0;
Mdl0 = arima('ARLags',1,'Constant',0,'Variance',VarMdl0);

% The structure of Mdl0 is the same as Mdl. However, every parameter is unknown, except for the restriction \alpha_1 = 0. These are equality constraints during estimation. You can interpret Mdl0 as an AR(1) model with the Gaussian innovations that have mean 0 and constant variance.
%
% Estimate the restricted model using the simulated data (y).

psI = 1:n;             % Presample indices
esI = (n + 1):(T + n); % Estimation sample indices

[EstMdl0,EstParamCov] = estimate(Mdl0,y(esI),...
    'Y0',y(psI),'E0',ep(psI),'V0',v(psI),'display','off');
phi10 = EstMdl0.AR{1};
alpha00 = EstMdl0.Variance.Constant;

% EstMdl0 contains the parameter estimates of the restricted model.
%
% lmtest requires the unrestricted model score evaluated at the restricted model estimates. The unrestricted model loglikelihood function is
% l(\phi_1,\alpha_0,\alpha_1) = \sum_{t=2}^T\left(-0.5\log(2\pi) - 0.5\log h_t - \frac{\varepsilon_t^2}{2h_t}\right),
% where \varepsilon_t = y_t - \phi_1y_{t-1}. The unrestricted gradient is
% \frac{\partial l (\phi_1,\alpha_0,\alpha_1)}{\partial \alpha} =
% \sum_{t=2}^T\frac{1}{2h_t}z_tf_t,
% where z_t = [1, \varepsilon_{t-1}^2] and f_t = \frac{\varepsilon_t^2}{h_t} - 1. The information matrix is
% I = \frac{1}{2h_t^2} \sum_{t=2}^T z_t'z_t.
% Under the null, restricted model, h_t = h_0 = \hat\alpha_0 for all t, where \hat\alpha_0 is the estimate from the restricted model analysis.
% Evaluate the gradient and information matrix under the restricted model. Estimate the parameter covariance by inverting the information matrix.

e = y - phi10*lagmatrix(y,1);
eLag1Sq = lagmatrix(e,1).^2;
h0 = alpha00;
ft = (e(esI).^2/h0 - 1);
zt = [ones(T,1),eLag1Sq(esI)]';

score0 = 1/(2*h0)*zt*ft;        % Score function
InfoMat0 = (1/(2*h0^2))*(zt*zt');
EstParamCov0 = inv(InfoMat0);   % Estimated parameter covariance
dof = 1;                        % Number of model restrictions

%
% Test the null hypothesis that \alpha_1=0 at the 5% significance level using lmtest.

[h,pValue] = lmtest(score0,EstParamCov0,dof)

% pValue is close to 0, which suggests that there is evidence to reject the restricted AR(1) model in favor of the unrestricted AR(1) model with an ARCH(1) variance.
% Copyright 2015 The MathWorks, Inc.