mfe_loadf_predb.m

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

function F=mfe_loadf_predB(L,stochcf,model_size,cal_beg_i,cal_end_i,for_beg_i,for_end_i,cal_beg_wday,for_beg_wday,X,VYEAR)
%MFE_LOADF_PREDB Auxiliary routine for MFE_LOADF.
%	F=MFE_LOADF_PREDB(...) returns a vector of Model B day-ahead 
%   forecasts of the daily time series L.
%
%   Input data:
%       L - daily load vector,
%       STOCHCF - stochastic component for the full period, i.e. for
%           (CAL_BEG_I:FOR_END_I),
%       MODEL_SIZE - ARMA(X) model size; type "help armax" at the command
%           line for format description,
%       CAL_BEG_I - index of the first value of the calibration period,
%       CAL_END_I - index of the last value of the calibration period,
%       FOR_BEG_I - index of the first value of the forecast period,
%       FOR_END_I - index of the last value of the forecast period,
%       CAL_BEG_WDAY - weekday of CAL_BEG_I,
%       FOR_BEG_WDAY - weekday of FOR_BEG_I,
%       X - exogenous variable, 
%       VYEAR - annual seasonal component for the full period.
%
%   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

% Compute 7-day seasonality (MA technique) in the calibration period
[L7,s7] = remst(L(cal_beg_i:cal_end_i),7,-1); 

% Remove 7-day seasonality (MA technique) in the full sample using 
% precomputed weekly component 's7'
L7 = remst(L(cal_beg_i:for_end_i),7,s7); 

% Define calibration and forecast periods
cal_period = 1+cal_end_i-cal_beg_i;
for_period = 1+for_end_i-for_beg_i;

% Compute day-ahead forecasts
y = zeros(for_period,1);
if size(model_size,2)==2
    disp(['ARMA [',num2str(model_size),']'])
    disp(['Processing ' num2str(for_period) ' data points ...'])
    for t = 1:for_period,
        % calibrate ARMA model to the stochastic component  
        M66 = armax(stochcf(t:t+cal_period-2),model_size); 
        % make a day-ahead prediction
        y66 = predict(stochcf(t:t+cal_period-1),M66);  
        y(t) = y66(end);
    end
else
    disp(['ARMAX [',num2str(model_size),']'])
    disp(['Processing ' num2str(for_period) ' data points ...'])
    for t = 1:for_period,
        % calibrate ARMAX model to the stochastic component
        M66 = armax([stochcf(t:t+cal_period-2) X(t:t+cal_period-2)],model_size); 
        % make a day-ahead prediction
        y66 = predict([stochcf(t:t+cal_period-1) X(t:t+cal_period-1)],M66);  
        y(t) = y66(end);
    end
end

% Shift the 7-day cycle if days of the week do not match
% convert to standard week (Sun, Mon, ..., Sat) = (1, 2, ..., 7)
s7wday = [s7((7-cal_beg_wday+2):7); s7(1:(7-cal_beg_wday+1))];
% shift 
s7wday = [s7wday(for_beg_wday:7); s7wday(1:(for_beg_wday-1))];

% invert the forecast vector
tt = (cal_end_i-cal_beg_i+1):(for_end_i-cal_beg_i);  
F = remst(L7(tt).*exp(y.*VYEAR(tt)),7,-s7wday);