startarx.m

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

function startarx(data,modelorder,preproc,ppc,Ndays,startd,endd,alphaci,obj1,obj2);
%STARTARX Auxiliary routine for MFE_PRICEF_WIN_ARX.
%	STARTARX(...) calculates and displays results of ARX model day-ahead 
%   forecasts.
%
%   Input data:
%       DATA - 5-column matrix (date, hour, price, load, load forecast),
%       MODELORDER - order of the ARX model,
%       PREPROC - spike preprocessing scheme (1 - limit, 2 - damped, 
%           3 - similar day, other - none),
%       PPC - spike preprocessing coefficient (spike is any value exceeding
%           data mean + PPS * data standard deviation),
%       NDAYS - number of days to be forecasted,
%       STARTD - index of the first day of the calibration period,
%       ENDD - index of the last day (in the first step) of the calibration 
%           period,
%       ALPHACI - confidence interval levels,
%       OBJ1 - pointer to confidence interval selection object,
%       OBJ2 - pointer to price cap selection object.
%
%   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

Naci = length(alphaci);
result = zeros(Ndays*24,8*Naci+5);

for j = 1:Ndays
    % display current day and number of days to be forecasted
    disp([j Ndays])
    
    % define price caps
    pricecap = 750;
    if data(endd+j*24,1)>=20000807,
        pricecap = 250;
    elseif data(endd+j*24,1)>=20000715,
        pricecap = 500;
    end
    
    % initialize 'result' matrix
    result((j-1)*24+1:j*24,1:3) = data(endd+(j-1)*24+1:endd+j*24,1:3);
    
    % perform one-day-ahead forecast
    resultaux = forecastarx(data(startd:endd+j*24,:),modelorder,preproc,ppc,alphaci);
    result((j-1)*24+1:j*24,4) = min(pricecap, resultaux(:,1));
    result((j-1)*24+1:j*24,5) = resultaux(:,1);
    
    % correct results for price cap
    for i=1:Naci
        result((j-1)*24+1:j*24,(i-1)*8+6) = min(pricecap, resultaux(:,(i-1)*4+2));
        result((j-1)*24+1:j*24,(i-1)*8+7) = min(pricecap, resultaux(:,(i-1)*4+3));
        result((j-1)*24+1:j*24,(i-1)*8+8) = min(pricecap, resultaux(:,(i-1)*4+4));
        result((j-1)*24+1:j*24,(i-1)*8+9) = min(pricecap, resultaux(:,(i-1)*4+5));
        result((j-1)*24+1:j*24,(i-1)*8+(10:13)) = resultaux(:,(i-1)*4+(2:5));
    end;
    
    % save temporary results
    save tempresult.txt result -ascii
end;

% print and plot results
[MDE,MWE,DRMSE,WRMSE,MeDE,MeWE] = ferrors(result(:,[3,4]));
mfe_pricef_errtables(MDE,MWE,DRMSE,WRMSE,MeDE,MeWE,data(endd+1:24:endd+Ndays*24,1));
mfe_pricef_plots(result,obj1,obj2);

% save final results
% format of 'result':
%   date, hour, actual price, price forecast, capped price forecast, 
%   [{lcl ucl gauss, lcl ucl npar}, 
%   price capped {lcl ucl gauss, lcl ucl npar}] as many times as there are
%   confidence levels (but no more than 3); a max of 29 columns in total.
save finalresult.txt result -ascii

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

function prediction=forecastarx(data,modelorder,preproc,ppc,alphaci);
%FORECASTARX Internally used by STARTARX.
%   PREDICTION=FORECASTARX(...) returns one-day-ahead point and interval 
%   forecasts. PREDICTION is a matrix with 24 rows. Point forecasts are in
%   the first column, interval forecasts in the following columns.
%
%   Input data:
%       DATA - 5-column matrix (date, hour, price, load, load forecast),
%       MODELORDER - order of the ARX model,
%       PREPROC - spike preprocessing scheme (1 - limit, 2 - damped, 
%           3 - similar day, other - none),
%       PPC - spike preprocessing coefficient (spike is any value exceeding
%           data mean + PPS * data standard deviation),
%       ALPHACI - confidence interval levels.

% weekday of starting point
datastart = data(1,1);
dayd = rem(datastart,100);
monthd = floor(rem(datastart,10000)/100);
yeard = floor(datastart/10000);
i = weekday(datestr(datenum(yeard,monthd,dayd,0,0,0),1))-1;

alphaci = (1+alphaci/100)/2;
Naci = length(alphaci);

N = length(data);
data = [data zeros(N,3)];
for j=1:24:N
    switch mod(i,7)
        case 6
            data(j:j+23,6)=1;   % Saturday
        case 0
            data(j:j+23,7)=1;   % Sunday
        case 1
            data(j:j+23,8)=1;   % Monday
    end;
    i=i+1;
end;

% spike preprocessing
if ismember(preproc,1:3)
    datamod=data(1:end-24,3);
    m1=mean(datamod);
    std1=std(datamod);
    limpp=m1+ppc*std1;
    indpp=find(datamod>limpp);
    switch preproc
        case 1 % limit
            datamod(indpp)=limpp;
        case 2 % damped
            datamod(indpp)=limpp*(1+log10(datamod(indpp)/limpp));
        case 3 % similar day
            dataaux=data(1:end-24,6:8);
            for nrind=indpp'
                if all(dataaux(nrind,:))==0
                    if nrind>24
                        datamod(nrind)=datamod(nrind-24,1);
                    else
                        datamod(nrind)=limpp*(1+log10(datamod(nrind,1)/limpp));
                    end
                else
                    if nrind>168
                        datamod(nrind)=datamod(nrind-168,1);
                    else
                        datamod(nrind)=limpp*(1+log10(datamod(nrind,1)/limpp));
                    end;
                end
            end;
    end;
    data(1:end-24,3)=datamod;
end;

% initialize 'prediction' matrix
prediction = zeros(24,4*Naci+1);

% compute forecasts
for hour = 1:24
    % preliminaries
    price = data(hour+168:24:end-24,3);
    price = log(price);
    mc = mean(price);
    price = price-mc;

    price_168 = data(hour:24:end-168,3);
    price_168 = log(price_168);
    price_168 = price_168-mean(price_168);

    price_min = data(169-24:end-24,3);
    price_min = reshape(price_min,24,length(price_min)/24);
    price_min = log(price_min);
    price_min = min(price_min)';
    price_min = price_min-mean(price_min);

    loadr = data(hour+168:24:end-24,4);
    loadr = log(loadr);
    loadr = loadr-median(loadr);

    loadf = data(hour+168:24:end,5);
    loadf = log(loadf);
    loadf = loadf-mean(loadf);
    
    % calibrate AR/ARX model
    MODEL = armax([price price_168(1:end-1) price_min(1:end-1) [loadr(1:end-1);loadf(end-1)] data(hour+168:24:end-24,6:8)],modelorder);
    data_pred = [[price;0],price_168,price_min,[loadr;loadf(end)],data(hour+168:24:end,6:8)];
    
    % perform forecast
    prog = predict(MODEL,data_pred);
    prediction(hour,1) = exp(prog(end)+mc);
    Npred = length(prog)-1;
    auxpred = sort(data_pred(1:Npred,1)-prog(1:Npred));
    
    % compute interval forecasts
    for i=1:Naci
        %gaussian interval
        prediction(hour,(i-1)*4+2)=prediction(hour,1)/exp(norminv(alphaci(i))*std(auxpred));
        prediction(hour,(i-1)*4+3)=prediction(hour,1)*exp(norminv(alphaci(i))*std(auxpred));
        %non-parametric interval
        prediction(hour,(i-1)*4+4)=prediction(hour,1)*exp(auxpred(ceil(Npred*(1-alphaci(i)))));
        prediction(hour,(i-1)*4+5)=prediction(hour,1)*exp(auxpred(floor(Npred*alphaci(i))));
    end;
end;