%% ARMA Model Selection
%
%% Clear the environment
% It is a good idea to clear all variables and close all figures before proceeding

clear all
close all
clc

%
%% Load the data
% The next step is to load the data. The next two lines set the maximum AR and MA lags.

load term
MAX_AR = 4;
MAX_MA = 4;

%
%% Estimation
% Model selection is challenging to automate. Instead, this code estimates all ARMAR models with p=0,1,...,4 and q=0,1,...4, and then looks at log likelihoods, AICs and BICs.

% Use a cell array to hold results from the models.  Cell arrays allow
% arbitrary data to be stored.
models = cell(MAX_AR+1,MAX_MA+1);

% Options to suppress display.  ARMA models are estimated using lsqnonlin
options =  optimset('lsqnonlin');
options.MaxIter = 1000;
options.Display='none';

% Setup arrays to hold values
LLs = zeros(MAX_AR,MAX_MA);
AICs = zeros(MAX_AR,MAX_MA);
BICs = zeros(MAX_AR,MAX_MA);

% Get the effective T since MAX_AR will be held back
T = length(term) - MAX_AR;
% Loop over AR and MA order
for i=0:MAX_AR
    for j=0:MAX_MA
        disp(['AR: ' num2str(i) ' MA: ' num2str(j)])
        % The MAX_AR term below enforces the holdbac, which is requires
        % when comparing AR models with different lag lengths.  Failing to
        % use this will produce log-likelihods base don different number of
        % observations.
        [p,ll,ser] = armaxfilter(term,1,1:i,1:j,[],[],options,MAX_AR);
        % Use a struct to store results
        model = struct();
        model.parameters = p;
        model.LL = ll;
        model.SER = ser;
        models{i+1,j+1} = model;
        % Store other values for use later
        LLs(i+1,j+1) = ll;
        AICs(i+1,j+1) = ll - 2 * (1+i+j);
        BICs(i+1,j+1) = ll - log(T) * (1+i+j);
    end
end

%
%% Plot the Log-likelihoods
% This produces a surface plot of the log-likelihoods. The AR=0 case has been eliminated since the log-lieklihoods of this model are very bad. Unsuprisingly, the largest model has the largest LL, although the gains appear to be small once the AR and MA order are both >= 1

% New figure
figure()
% Use surf to produce a surface plot
surf(0:4,1:4,LLs(2:5,:))
% Labels - X is the first (columns of data), Y is the second (rows)
xlabel('MA Order')
ylabel('AR Order')
% Title
title('Log-likelihoods')
disp('LL difference from highest')
disp(LLs(2:5,:) - max(max(LLs)))

%
%% Plot the AICs
% The AICs are plotted in the same matter. The maximum value is also computed.

% New figure
figure()
% Use surf to produce a surface plot
surf(0:4,1:4,AICs(2:5,:))
% Labels - X is the first (columns of data), Y is the second (rows)
xlabel('MA Order')
ylabel('AR Order')
% Title
title('AIC')
% Find the model which has the maximum value
[AIC_AR,AIC_MA] = find(AICs==max(max(AICs)));
disp('AIC Maximum at:')
disp([AIC_AR - 1 AIC_MA - 1]);

%
%% Plot the BICs
% The BICs are plotted in the same matter. The maximum value is also computed.

% New figure
figure()
% Use surf to produce a surface plot
surf(0:4,1:4,BICs(2:5,:))
% Labels - X is the first (columns of data), Y is the second (rows)
xlabel('MA Order')
ylabel('AR Order')
% Title
title('BIC')
[BIC_AR,BIC_MA] = find(BICs==max(max(BICs)));
disp('BIC Maximum at:')
disp([BIC_AR - 1 BIC_MA - 1]);