mfe_distrib_aux.m

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

function mfe_distrib_aux(x,dists,edft,radiobut,rems);
%MFE_DISTRIB_AUX Auxiliary routine for MFE_DISTRIB.
%	MFE_DISTRIB_AUX(X,DISTS,EDFT,RADIOBUT,REMS) runs the computations and
%	displays results for Case Study 2.6.3 in [1]. 
%
%   Input data:
%       X - data vector,
%       DISTS - handle to the distributions checkboxes,
%       EDFT - handle to the edf checkbox,
%       RADIOBUT - handle to the preprocessing radiobuttons,
%       REMS - handle to the 'remove weekly cycle' checkbox.
%
%   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


disp('---------------------------------------------------------------------------------------------------')
disp('                                 Parameters                 ||                Test values          ')
disp('---------------------------------------------------------------------------------------------------')
disp('               alpha      sigma, delta      beta      mu    || Anderson-Darling      Kolmogorov    ')
disp('---------------------------------------------------------------------------------------------------')

% Remove weekly component
if get(rems,'value')
    x = remst(x,7,-2);
end;

% Plot empirical data
h=figure(5);
set(h,'name','Data','numbertitle','off');
switch radiobut
    case 'dif'
        x = diff(x);
        plot(x,'b');
        ylabel('Price changes');
    case 'ret'
        x = 100*logret(x);
        plot(x,'b');
        ylabel('Returns [%]');
    otherwise
        plot(x,'b');
        ylabel('Prices');
end
xlabel('Days')

% Compute the empirical cdf
[X,Y] = empcdf(x);
% remove the infimum of the support
X = X(2:end); Y = Y(2:end);

% Compute and plot edf
h=figure(2);
set(h,'name','CDF','numbertitle','off');
plot(X,Y,'b.')
xlabel('x')
ylabel('CDF(x)')
hold on
xplus=find(X>0);
lxplus=sum(xplus);
if lxplus>0
    h=figure(3);
    set(h,'name','Right tail','numbertitle','off');
    loglog(X(xplus),1-Y(xplus),'b.')
    set(gca,'ylim',[min(1-Y(xplus(1:end-1)))/2 max(1-Y(xplus))*2])
    xlabel('x')
    ylabel('1-CDF(x)')
    hold on
end;
xminus=find(X<0);
lxminus=sum(xminus);
if lxminus>0
    h=figure(4);
    set(h,'name','Left tail','numbertitle','off');
    loglog(-X(xminus),Y(xminus),'b.')
    set(gca,'ylim',[min(Y(xminus))/2 max(Y(xminus))*2])
    xlabel('-x')
    ylabel('CDF(x)')
    hold on
end;
ld={'EDF'};

% Compute and plot Gaussian fit
if get(dists(1),'value')
    params=[mean(x),std(x)];
    ts=[];
    if get(edft,'value')
        [A2,K]=edftests(x,params,'Gaussian');
    end;
    rprint(params,[A2,K],'Gaussian    ')   
    figure(2)
    plot(X,normcdf(X,params(1),params(2)),'c')
    if lxplus>0
        figure(3)
        loglog(X(xplus),1-normcdf(X(xplus),params(1),params(2)),'c')
    end
    if lxminus>0
        figure(4)
        loglog(-X(xminus),normcdf(X(xminus),params(1),params(2)),'c')
    end;
    ld={ld{:},'Gaussian'};
end;

% Compute and plot hyperbolic fit
if get(dists(2),'value')
    params=hypest(x);
    ts=[];
    if get(edft,'value')
        [A2,K]=edftests(x,params,'hyperbolic');
    end;
    rprint([params(1) params(3) params(2) params(4)],[A2,K],'Hyperbolic  ') 
    figure(2)
    plot(X,hypcdf(X,params(1),params(2),params(3),params(4)),'r')
    if lxplus>0
        figure(3)
        plot(X(xplus),1-hypcdf(X(xplus),params(1),params(2),params(3),params(4)),'r')
    end
    if lxminus>0
        figure(4)
        plot(-X(xminus),hypcdf(X(xminus),params(1),params(2),params(3),params(4)),'r')
    end
    ld={ld{:},'Hyperbolic'};
end;

% Compute and plot NIG fit
if get(dists(3),'value')
    params=nigest(x);
    ts=[];
    if get(edft,'value')
        [A2,K]=edftests(x,params,'NIG');
    end;
    rprint([params(1) params(3) params(2) params(4)],[A2,K],'NIG         ')
    figure(2)
    plot(X,nigcdf(X,params(1),params(2),params(3),params(4)),'m')
    if lxplus>0
        figure(3)
        loglog(X(xplus),1-nigcdf(X(xplus),params(1),params(2),params(3),params(4)),'m')
    end
    if lxminus>0
        figure(4)
        loglog(-X(xminus),nigcdf(X(xminus),params(1),params(2),params(3),params(4)),'m')
    end;
    ld={ld{:},'NIG'};
end;

% Compute and plot alpha-stable fit
if get(dists(4),'value')
    % use regression estimator of Koutrouvelis (1980)
    [alpha,sigma,beta,mu]=stabreg(x);
    params=[alpha,sigma,beta,mu];
    ts=[];
    if get(edft,'value')
        [A2,K]=edftests(x,params,'stable');
    end;
    rprint(params,[A2,K],'Alpha-stable')
    figure(2)
    plot(X,stabcdf(X,alpha,sigma,beta,mu),'k')
    if lxplus>0
        figure(3)
        loglog(X(xplus),1-stabcdf(X(xplus),alpha,sigma,beta,mu),'k')
    end
    if lxminus>0
        figure(4)
        loglog(-X(xminus),stabcdf(X(xminus),alpha,sigma,beta,mu),'k')
    end;
    ld={ld{:},'Alpha-stable'};
end;

% Add legends to the figures
figure(2)
legend(ld,4);
set(gca,'ylim',[0,1]);
hold off
if lxplus>0
    figure(3)
    legend(ld,3);
    hold off
end;
if lxminus>0
    figure(4)
    legend(ld,3);
    hold off
end;

%=========================================================================
% Internally used routine(s)
%=========================================================================

function rprint(params,edft,dist);
%RPRINT Internally used by MFE_DISTRIB_AUX.
%   RPRINT(PARAMS,EDFT,DIST) Displays the parameter estimates and edf
%   statistics.
%
%   Input data:
%       PARAMS - parameter estimates vector
%       EDFT - edf statistics vector
%       DIST - distribution name (string)

switch dist
    case 'Gaussian    '
        params=['              ',sprintf('%7.4f',params(2)),'                ',sprintf('%7.4f',params(1))];       
otherwise
        params=[sprintf('%7.4f',params(1)),'       ',sprintf('%7.4f',params(2)),'      ',sprintf('%7.4f',params(3)),'   ',sprintf('%7.4f',params(4))];
end;
if ~isempty(edft)
    edft=[sprintf('%7.4f',edft(1)),'             ',sprintf('%7.4f',edft(2))];
end;
disp([dist,'  ',params,'  ||    ',edft])