mfe_loadf_crit.m

Modeling and Forecasting Electricity Loads and Prices: A Statistical Approach" by Rafa&#322 Weron, p

function mfe_loadf_crit(stochc,maxp,maxq,crit);
%MFE_LOADF_CRIT Auxiliary routine for MFE_LOADF.
%	MFE_LOADF_CRIT(STOCHC,MAXP,MAXQ,CRIT) searches for the best ARMA(p,q)
%	model (with p=1,...,MAXP and q=1,...,MAXQ) of the stochastic component
%	STOCHC in terms of the information criterion CRIT. CRIT=3 stands for
%	AICC, while CRIT=4 for BIC. Results are displayed in the command
%	window in two tables:
%   1. Best ARMA(p,q) model for a given AR order p=1,...,MAXP
%   2. Best 5 models
%
%   Reference(s):
%   [1] R.Weron (2007) "Modeling and Forecasting Electricity Loads and 
%   Prices: A Statistical Approach", Wiley, Chichester.

%   Written by Adam Misiorek and Rafal Weron (2006.09.22)
%   Copyright (c) 2006 by Rafal Weron

% Remove the mean
stochc = stochc-mean(stochc);
n = length(stochc);

rtab = [];
for i = 0:maxp
    for j = 0:maxq
        if i+j>0
            model = armax(stochc,[i,j]);
            Lpred = predict(model,stochc);
            res = stochc-Lpred;
            s2 = var(res);
            ACF = armaacvf(model.a,model.c,n)';
            V = zeros(n);
            for k = 1:n
                V(k,:) = ACF(abs((1:n)-k)+1);
            end;
            % compute -2*log-likelihood
            % 'abs' in 'abs(det(V)' is in case of numerical inaccuracies
            m2logL = n*log(2*pi*s2)+log(abs(det(V)))+stochc'*inv(V)*stochc/s2; 
        else
            res = stochc;
            s2 = var(res);
            m2logL = n*(log(2*pi*s2)+1);
        end;
        % [AR order, MA order, AICC, BIC]
        rtab = [rtab; i j m2logL+2*n*(i+j+1)/(n-i-j-2) m2logL+(i+j+1)*log(n)];
    end;
end;

switch crit
    case 3,
        disp('-----------------------------------------------------------------------')
        disp('        AR order      MA order      AICC*          BIC')
        disp('-----------------------------------------------------------------------')
    case 4,
        disp('-----------------------------------------------------------------------')
        disp('        AR order      MA order      AICC          BIC*')
        disp('-----------------------------------------------------------------------')
end        

disp('-----------------------------------------------------------------------')
disp('                    Best model for each AR order ')
disp('-----------------------------------------------------------------------')

for i=0:maxp
    pom=sortrows(rtab(i*(maxq+1)+1:(i+1)*(maxq+1),:),crit);
    disp(sprintf('            %d             %d       %8.2f      %8.2f',pom(1,:)))
end;

disp('-----------------------------------------------------------------------')
disp('                           Best 5 models ')
disp('-----------------------------------------------------------------------')

pom = sortrows(rtab,crit);
l = size(pom);
for i=1:min(5,l(1))
    disp(sprintf('            %d             %d       %8.2f      %8.2f',pom(i,:)))
end;