startgarch.m

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

function startgarch(data,modelorder,Ndays,startd,endd,alphaci,obj1,obj2,alg);
%STARTGARCH Auxiliary routine for MFE_PRICEF_WIN_ARXG.
%	STARTARX(...) calculates and displays results of ARX-GARCH (or 
%   ARX+GARCH) model day-ahead forecasts. In the ARX-GARCH model all
%   parameters are estimated simultanously. In the ARX+GARCH model first an
%   ARX model is calibrated, then the residuals are fitted by a GARCH
%   model.
%
%   Input data:
%       DATA - 5-column matrix (date, hour, price, load, load forecast),
%       MODELORDER - order of the ARX model,
%       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,
%       ALG - model type switch: ALG=1 -> ARX-GARCH, else -> ARX+GARCH.
%
%   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 for the selected model
    if alg==1, %'ARX-GARCH',
        resultaux=forecastgarch(data(startd:endd+j*24,:),modelorder,alphaci);
    else    %'ARX + GARCH'
        resultaux=forecastarxgarch(data(startd:endd+j*24,:),modelorder,alphaci);
    end
    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=forecastgarch(data,modelorder,alphaci);
%FORECASTGARCH Internally used by STARTGARCH.
%   PREDICTION=FORECASTGARCH(...) returns one-day-ahead point and interval 
%   forecasts of the ARX-GARCH model. 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,
%       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;

% initialize 'prediction' matrix
prediction = zeros(24,5);

% 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);
    
    if modelorder(4) % ARX-GARCH
        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 ARX-GARCH model
        [Coeff,Errors,LLF,Innovations,Sigma,Summary] = garchfit(garchset('display','off','r',modelorder(1),'p',modelorder(2),'q',modelorder(3)),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)]);
        % perform forecast
        [aux,prog] = garchpred(Coeff,price,1,[price_168 price_min [loadr;loadf(end)] data(hour+168:24:end,6:8)]);
    else % AR-GARCH
        % calibrate AR-GARCH model
        [Coeff,Errors,LLF,Innovations,Sigma,Summary] = garchfit(garchset('display','off','r',modelorder(1),'p',modelorder(2),'q',modelorder(3)),price,[price_168(1:end-1) price_min(1:end-1) data(hour+168:24:end-24,6:8)]);
        % perform forecast
        [aux,prog] = garchpred(Coeff,price,1,[price_168 price_min data(hour+168:24:end,6:8)]);
    end;    
    prediction(hour,1) = exp(prog+mc);
    Npred = length(Innovations);
    Innovations = sort(Innovations);

    % compute interval forecasts
    for i=1:Naci
        %gaussian interval
        prediction(hour,(i-1)*4+2)=prediction(hour,1)/exp(norminv(alphaci(i))*std(Innovations));
        prediction(hour,(i-1)*4+3)=prediction(hour,1)*exp(norminv(alphaci(i))*std(Innovations));
        %non-parametric interval
        prediction(hour,(i-1)*4+4)=prediction(hour,1)*exp(Innovations(ceil(Npred*(1-alphaci(i)))));
        prediction(hour,(i-1)*4+5)=prediction(hour,1)*exp(Innovations(floor(Npred*alphaci(i))));
    end;
end;


function prediction=forecastarxgarch(data,modelorder,alphaci);
%FORECASTARXGARCH Internally used by STARTGARCH.
%   PREDICTION=FORECASTARXGARCH(...) returns one-day-ahead point and interval 
%   forecasts of the ARX+GARCH model. 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,
%       ALPHACI - confidence interval levels.

% weekday of start 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;

% initialize 'prediction' matrix
prediction=zeros(24,5);

% 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(1),1,1,modelorder(4),1,1,1,0,0,0,0,0,0,0]);
    data_pred = [[price;0],price_168,price_min,[loadr;loadf(end)],data(hour+168:24:end,6:8)];

    % perform forecast
    proga = predict(MODEL,data_pred);

    Npred=length(price);
    resg=price-proga(1:Npred);

    % calibrate GARCH model to residuals
    [Coeff,Errors,LLF,Innovations,Sigma,Summary]=garchfit(garchset('display','off','p',modelorder(2),'q',modelorder(3)),resg);

    % perform forecast
    [aux,progg]=garchpred(Coeff,resg,1);

    % combine AR/ARX and GARCH forecasts
    prediction(hour,1)=exp(proga(end)+progg+mc);
    Innovations=sort(Innovations);

    % compute interval forecasts
    for i=1:Naci
        %gaussian interval
        prediction(hour,(i-1)*4+2)=prediction(hour,1)/exp(norminv(alphaci(i))*std(Innovations));
        prediction(hour,(i-1)*4+3)=prediction(hour,1)*exp(norminv(alphaci(i))*std(Innovations));
        %non-parametric interval
        prediction(hour,(i-1)*4+4)=prediction(hour,1)*exp(Innovations(ceil(Npred*(1-alphaci(i)))));
        prediction(hour,(i-1)*4+5)=prediction(hour,1)*exp(Innovations(floor(Npred*alphaci(i))));
    end;
end;