%% Test Among Multiple Nested Model Specifications
% Assess model specifications by testing down among multiple restricted models using simulated data. The true model is the ARMA(2,1)
% y_t = 3 + 0.9y_{t-1}-0.5y_{t-2}+\varepsilon_t+0.7\varepsilon_{t-1},
% where \varepsilon_t is Gaussian with mean 0 and variance 1.
% Specify the true ARMA(2,1) model, and simulate 100 response values.

TrueMdl = arima('AR',{0.9,-0.5},'MA',0.7,...
    'Constant',3,'Variance',1);
T = 100;
rng(1); % For reproducibility 
y = simulate(TrueMdl,T);

%
% Specify the unrestricted model and the names of the candidate models for testing down.

UMdl = arima(2,0,2);
RMdlNames = {'ARMA(2,1)','AR(2)','ARMA(1,2)','ARMA(1,1)',...
    'AR(1)','MA(2)','MA(1)'};

% UMdl is the unrestricted, ARMA(2,2) model. RMdlNames is a cell array of strings containing the names of the restricted models.
%
% Fit the unrestricted model to the simulated data.

[UEstMdl,UEstParamCov] = estimate(UMdl,y,'Display','off');

% UEstMdl is the fitted, unrestricted model, and UEstParamCov is the estimated parameter covariance matrix.
% The unrestricted model has six parameters. To construct the restriction function and its Jacobian, you must know the order of the parameters in UEstParamCov. For this arima model, the order is [c,\phi_1,\phi_2,\theta_1\theta_2,\sigma^2].
%
% Each candidate model corresponds to a restriction function. Put the restriction function vectors into separate cells of a cell vector.

rf1 = UEstMdl.MA{2};                          % ARMA(2,1)
rf2 = cell2mat(UEstMdl.MA)';                  % AR(2)
rf3 = UEstMdl.AR{2};                          % ARMA(1,2)
rf4 = [UEstMdl.AR{2};UEstMdl.MA{2}]';         % ARMA(1,1)
rf5 = [UEstMdl.AR{2};cell2mat(UEstMdl.MA)'];  % AR(1)
rf6 = cell2mat(UEstMdl.AR)';                  % MA(2)
rf7 = [cell2mat(UEstMdl.AR)';UEstMdl.MA{2}];  % MA(1)
r = {rf1;rf2;rf3;rf4;rf5;rf6;rf7};

% r is a 7-by-1 cell vector of vectors corresponding to the restriction function for the candidate models.
%
% Put the Jacobian of each restriction function into separate, corresponding cells of a cell vector. The order of the elements in the Jacobian must correspond to the order of the elements in UEstParamCov.

J1 = [0 0 0 0 1 0];                           % ARMA(2,1)   
J2 = [0 0 0 1 0 0; 0 0 0 0 1 0];              % AR(2)      
J3 = [0 1 0 0 0 0];                           % ARMA(1,2)  
J4 = [0 1 0 0 0 0; 0 0 0 0 1 0];              % ARMA(1,1)  
J5 = [0 1 0 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0]; % AR(1)      
J6 = [1 0 0 0 0 0; 0 1 0 0 0 0];              % MA(2)      
J7 = [1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 0 0 1 0]; % MA(1)      
R = {J1;J2;J3;J4;J5;J6;J7};

% R is a 7-by-1 cell vector of vectors corresponding to the restriction function for the candidate models.
%
% Put the estimated parameter covariance matrix in each cell of a 7-by-1 cell vector.

EstCov = cell(7,1); % Preallocate
for j = 1:length(EstCov)
    EstCov{j} = UEstParamCov;
end

%
% Apply the Wald test at a 1% significance level to find the appropriate, restricted model specifications.

alpha = .01;
h = waldtest(r,R,EstCov,alpha);
RestrictedModels = RMdlNames(~h)

% RestrictedModels lists the most appropriate restricted models.
% You can test down again, but use ARMA(2,1) as the unrestricted model. In this case, you must remove MA(2) from the possible restricted models.
% Copyright 2015 The MathWorks, Inc.