starttarx.m

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

function starttarx(data,aro,ext,thrvar,Ndays,startd,endd,alphaci,obj1,obj2);
%STARTTARX Auxiliary routine for MFE_PRICEF_WIN_TARX.
%	STARTTARX(...) calculates and displays results of TARX model day-ahead 
%   forecasts.
%
%   Input data:
%       DATA - 5-column matrix (date, hour, price, load, load forecast),
%       ARO - AR order,
%       EXT - exogenous variable switch: EXT=0 -> TAR, else -> TARX,
%       THRVAR - threshold variable switch: THRVAR=1 -> load, else -> price,
%       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=forecasttarx(data(startd:endd+j*24,:),aro,ext,thrvar,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=forecasttarx(data,aro,ext,thrvar,alphaci);
%FORECASTTARX Internally used by STARTTARX.
%   PREDICTION=FORECASTTARX(...) 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),
%       ARO - AR order,
%       EXT - exogenous variable switch: EXT=0 -> TAR, else -> TARX,
%       THRVAR - threshold variable switch: THRVAR=1 -> load, else -> price,
%       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-24-168+hour,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);
    
    % threshold variable
    if thrvar % select load as threshold variable
        weekdiff=data(1:end-24,4);
    else % select price as threshold variable
        weekdiff=data(1:end-24,3);
    end;
    weekdiff=mean(reshape(weekdiff,24,length(weekdiff)/24));
    weekdiff=weekdiff(8:end)-weekdiff(1:end-7);
    weekdiff=weekdiff'-mean(weekdiff);
    
    if ext % TARX
        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 TARX model
        [MODEL,resid] = modeltarx(price(2:end),weekdiff(1:end-1),[price_168(2:end-1) price_min(2:end-1) [loadr(2:end-1);loadf(end-1)] data(hour+24+168:24:end-24,6:8)],1:aro,1:aro);
        % perform forecast
        prediction(hour,1) = exp(tarxpred(price,weekdiff(end),[price_168(end) price_min(end) loadf(end) data(end-24+hour,6:8)],MODEL,1:2,1:2,1)+mc);
    else % TAR
        % calibrate TAR model
        [MODEL,resid] = modeltarx(price(2:end),weekdiff(1:end-1),[price_168(2:end-1) price_min(2:end-1) data(hour+24+168:24:end-24,6:8)],1:aro,1:aro);
        % perform forecast
        prediction(hour,1) = exp(tarxpred(price,weekdiff(end),[price_168(end) price_min(end) data(end-24+hour,6:8)],MODEL,1:2,1:2,1)+mc);
    end;    
    auxpred = sort(resid);
    Npred = length(auxpred);

    % 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;


function [estpar,resid]=modeltarx(x,z,outv,p1,p2);
%MODELTARX Internally used by FORECASTTARX.
%   [ESTPAR,RESID]=MODELTARX(...) returns parameter estimates and residuals 
%   of a 2-state TAR/TARX model.
%
%   Input data:
%       X - historical values of the forecasted process,
%       Z - historical values of the threshold variable,
%       OUTV - historical values of the exogenous variable(s),
%       P1 - lags of the first regime AR process, e.g. P1=1:ARO,
%       P2 - lags of the second regime AR process.

% fix threshold level of the threshold variable
qint = 0; 

% define sizes of AR processes
lp1 = length(p1);
lp2 = length(p2);
start = max(max(p1),max(p2));

% define regimes
zz = z(start+1:end)<=qint;
ind1 = find(zz);
ind2 = find(~zz);

% initialize parameter estimates vectors 
Noutvar = size(outv,2);
b1 = NaN*ones(lp1+Noutvar,1);
b2 = NaN*ones(lp2+Noutvar,1);

% return if too few values in either of the regimes
if size(ind1,1)<=lp1+Noutvar+1|size(ind2,1)<=lp2+Noutvar+1
    s2 = realmax;
    return;
end;

xx = x(start+1:end);
outvx = outv(start+1:end,:);

% select first regime data
data1 = [xx(ind1),x(start+ind1*ones(1,lp1)-ones(size(ind1))*p1) outvx(ind1,:)];

% calibrate first regime process
model1 = armax(data1,[0 ones(1,lp1+Noutvar) zeros(1,lp1+Noutvar+1)]);
b1 = model1.b;

% forecast first regime values
r1 = x(ind1)-predict(model1,data1);

% select second regime data
data2 = [xx(ind2),x(start+ind2*ones(1,lp2)-ones(size(ind2))*p2) outvx(ind2,:)];

% calibrate second regime process
model2 = armax(data2,[0 ones(1,lp2+Noutvar) zeros(1,lp2+Noutvar+1)]);
b2 = model2.b;

% forecast second regime values
r2 = x(ind2)-predict(model2,data2);

% compute residuals
resid = [r1;r2];

% define parameter estimates vector
estpar = [qint b1' b2'];


function y=tarxpred(x,z,outv,params,p1,p2,n);
%TARXPRED Internally used by FORECASTTARX.
%   Y=TARXPRED(...,N) returns N-step-ahead forecast of a 2-state TAR/TARX 
%   model.
%
%   Input data:
%       X - historical values of the forecasted process,
%       Z - future values of the threshold variable,
%       OUTV - values of the exogenous variable(s),
%       PARAMS - model parameters,
%       P1 - lags of the first regime AR process, e.g. P1=1:ARO,
%       P2 - lags of the second regime AR process,
%       N - number of forecast steps.

Noutvar = size(outv,2);
lp1 = length(p1);
lp2 = length(p2);
q = params(1);
a1 = params(2:lp1+1);
coeff1 = params(lp1+2:lp1+Noutvar+1);
a2 = params(lp1+Noutvar+2:lp1+Noutvar+lp2+1);
coeff2 = params(lp1+Noutvar+lp2+2:end);
T = length(x);
y = [x;zeros(n,1)];
for i=1:n
    if z(i)<=q % first regime
        y(T+i) = a1*y(T+i-p1)+coeff1*outv(i,:)';
    else % second regime
        y(T+i) = a2*y(T+i-p2)+coeff2*outv(i,:)';
    end;
end;
y = y(T+1:end);