%% Choose the Best AR Model Specification
% Compare AR model specifications for a simulated response series using  lmtest.
% Consider the AR(3) model:
% y_t = 1 + 0.9y_{t-1}-0.5y_{t-2}+0.4y_{t-3}+\varepsilon_t,
% where \varepsilon_t is Gaussian with mean 0 and variance 1. Specify this model using arima.

Mdl = arima('Constant',1,'Variance',1,'AR',{0.9,-0.5,0.4});

% Mdl is a fully specified, AR(3) model.
%
% Simulate presample and effective sample responses from Mdl.

T = 100;
rng(1);               % For reproducibility
n = max(Mdl.P,Mdl.Q); % Number of presample observations
y = simulate(Mdl,T + n);

% y is a a random path from Mdl that includes presample observations.
%
% Specify the restricted model:
% y_t = c + \phi_1y_{t-1}+\phi_2y_{t-2}+\varepsilon_t,
% where \varepsilon_t is Gaussian with mean 0 and variance \sigma^2.

Mdl0 = arima(3,0,0);
Mdl0.AR{3} = 0;

% The structure of Mdl0 is the same as Mdl. However, every parameter is unknown, except that \phi_3 = 0. This is an equality constraint during estimation.
%
% Estimate the restricted model using the simulated data (y).

[EstMdl0,EstParamCov] = estimate(Mdl0,y((n+1):end),...
    'Y0',y(1:n),'display','off');
phi10 = EstMdl0.AR{1};
phi20 = EstMdl0.AR{2};
phi30 = 0;
c0 = EstMdl0.Constant;
phi0 = [c0;phi10;phi20;phi30];
v0 = EstMdl0.Variance;

% EstMdl0 contains the parameter estimates of the restricted model.
%
% lmtest requires the unrestricted model score evaluated at the restricted model estimates. The unrestricted model gradient is
% \frac{\partial l (\phi_1,\phi_2,\phi_3,c,\sigma^2;y_t,...,y_{t-3})}{\partial c} =
% \frac{1}{\sigma^2}(y_t - c - \phi_1y_{t-1}-\phi_2y_{t-2}-\phi_3y_{t-3})
% \frac{\partial l (\phi_1,\phi_2,\phi_3,c,\sigma^2;y_t,...,y_{t-3})}{\partial \phi_j} =
% \frac{1}{\sigma^2}(y_t - c - \phi_1y_{t-1}-\phi_2y_{t-2}-\phi_3y_{t-3})y_{t-j}
% \frac{\partial l (\phi_1,\phi_2,\phi_3,c,\sigma^2;y_t,...,y_{t-3})}{\partial \sigma^2} =
% -\frac{1}{2\sigma^2}+\frac{1}{2\sigma^4}(y_t - c - \phi_1y_{t-1}-\phi_2y_{t-2}-\phi_3y_{t-3})^2.

MatY = lagmatrix(y,1:3);
LagY = MatY(all(~isnan(MatY),2),:);
cGrad = (y((n+1):end)-[ones(T,1),LagY]*phi0)/v0;
phi1Grad = ((y((n+1):end)-[ones(T,1),LagY]*phi0).*LagY(:,1))/v0;
phi2Grad = ((y((n+1):end)-[ones(T,1),LagY]*phi0).*LagY(:,2))/v0;
phi3Grad = ((y((n+1):end)-[ones(T,1),LagY]*phi0).*LagY(:,3))/v0;
vGrad = -1/(2*v0)+((y((n+1):end)-[ones(T,1),LagY]*phi0).^2)/(2*v0^2);
Grad = [cGrad,phi1Grad,phi2Grad,phi3Grad,vGrad]; % Gradient matrix

score = sum(Grad)'; % Score under the restricted model

%
% Evaluate the unrestricted parameter covariance estimator using the restricted MLEs and the outer product of gradients (OPG) method.

EstParamCov0 = inv(Grad'*Grad);
dof = 1; % Number of model restrictions

%
% Test the null hypothesis that \phi_3=0 at a 1% significance level using lmtest.

[h,pValue] = lmtest(score,EstParamCov0,dof,0.1)

% pValue is close to 0, which suggests that there is strong evidence to reject the restricted, AR(2) model in favor of the unrestricted, AR(3) model.
% Copyright 2015 The MathWorks, Inc.