Home - Lectures - Publications - Software - Students

 

MFE Toolbox ver. 1.0.1

by

Rafał Weron

with great help of Jakub Jurdziak and Adam Misiorek

Prepared for use with "Modeling and Forecasting Electricity Loads and Prices: A Statistical Approach" by Rafał Weron, published by John Wiley and Sons, 2006.

Copyright (c) 2006-2007 by Rafał Weron

Date of last revision: 2007-11-05


Contents

Part I: The Toolbox

1. Installation
2.
Using the MFE Toolbox functions
3.
Toolbox contents

Part II: Case Studies

4. Distributions (Case Study 2.6.3)
5.
Load forecasting (Case Studies 3.4.7 and 3.4.10)
6.
Price forecasting (Case Studies 4.3.7, 4.3.8 and 4.3.11)
7.
Regime switching models (Case Study 4.4.5)


1. Installation

The MFE Toolbox functions are written in MATLAB 7.0, and should therefore work in that or any later version of MATLAB. Most of the functions have been also tested in MATLAB 6.5. Note, that some of the results may differ from those presented in the book if a later version of MATLAB is used.

To install the software, you should copy the MFE Toolbox folder as a folder in your MATLAB 'works' folder, and then activate the path to this new folder so that MATLAB knows it is there. Most of the functions should work if you have the Garch and System Identification Toolboxes in additional to the basic MATLAB package; however, a small number draw on functions in other MATLAB Toolboxes (e.g., the Statistics Toolbox) and will only work if the relevant other MATLAB Toolbox is installed.

Back to Contents


2. Using the MFE Toolbox functions

MATLAB functions or scripts (i.e. series of MATLAB commands) can be run from the MATLAB command line, which is given by the prompt '>>' in the MATLAB Command Window, i.e.:

>>

Scripts can be executed by simply typing in the relevant file name in the command line. For example, if you wish to run Case Study 4.4.5 demo script, you can do so by typing in the following command:

>> mfe_regime

Functions, on the other hand, generally require input parameters. Care must be taken to ensure that the input data is in appropriate form. To check data format requirements type '>> help <function name>' at the command line. For instance, to see the help screen of the 'remst.m' function type in the following command:

>> help remst

In return MATLAB will display:

REMST Remove seasonal component and trend.
Y = REMST(X,D,DEGREE) returns time series Y with removed polynomial
trend of degree DEGREE and seasonal components of period D from vector
X. REMST uses the moving average technique (see [1], Section 2.4.3).
Y = REMST(X,D,DEGREE,1) glues the first and last D/2 values to the end
and the beginning of X, respectively (recommended for short time
series, i.e. of length 2-3 D).
[Y,S,P] = REMST(X,D,DEGREE) returns seasonal component S and
polynomial coefficients P.

If DEGREE==-1 then only seasonal components are removed.
If DEGREE==-2 then only seasonal components are removed and the minimum
values of the original and deseasonalized series are aligned.
If DEGREE is a vector it is treated as a one-period (i.e. of length D)
seasonal component and is subtracted from X.

Reference(s):
[1] R.Weron (2007) 'Modeling and Forecasting Electricity Loads and
Prices: A Statistical Approach', Wiley, Chichester.

Thus, the following commands (i) define a seasonal series of length 100 and period 4, (ii) decompose it using the moving average technique (see [1], Section 2.4.3), (iii) plot the original and deseasoanlized series and (iv) display the seasonal component in the command window:

>> x=sin((1:100)'*pi/2)+rand(100,1)/10;
>> [xdes,season]=remst(x,4,-1);
>> plot([x xdes]); legend('Original series x','Deseasonalized series xdes')
>> s'

ans =

1.0000 0.0000 -1.0000 -0.0000

Back to Contents


3. Toolbox contents

The MFE Toolbox v. 1.0 consists of 58 functions, scripts and data files. Many functions include internally used routines, hence, the total number of functions is much larger. The files are grouped into seven categories:

1. Time series,
2. Distributions,
3. Tests and goodness-of-fit functions,
4. Markov regime switching (MRS) models,
5. GUI functions,
6. Demos,
7. Data files.

The Time series group includes functions for time series decomposition and transformation:

  • armaacvf - Autocovariance function of an ARMA process.
  • average - Weighted average.
  • decompA - Differencing-smoothing daily data decomposition.
  • decompB - Moving average with rolling volatility daily data decomposition.
  • logret - Logarithmic returns.
  • mainfcoeff - Coefficients of a MA(inf) process corresponding to ARMA(p,q).
  • periodog - Periodogram of a time series.
  • remmed - Remove mean- or median-based seasonal component.
  • remsin - Remove annual sinusoidal component and trend from daily data.
  • remst - Remove seasonal component and trend.
  • rollingvol - Annual rolling volatility.
  • volaplot - Volatility plot.

The Distributions group mainly covers estimation and cdf/pdf formulas of three popular heavy-tailed distributions: hyperbolic, NIG and (alpha-)stable. The full list follows:

  • empcdf - Empirical cumulative distribution function (cdf).
  • hypcdf - Hyperbolic cumulative distribution function (cdf).
  • hypest - Estimate parameters of the hyperbolic distribution.
  • hyploglik - Hyperbolic log-likelihood function.
  • hyppdf - Hyperbolic probability density function (pdf).
  • nigcdf - NIG cumulative distribution function (cdf).
  • nigest - Estimate parameters of the NIG distribution.
  • nigloglik - NIG log-likelihood function.
  • nigpdf - NIG probability density function (pdf).
  • stabcdf - (Alpha-)stable cumulative distribution function (cdf).
  • stabcull - Quantile parameter estimates of a stable distribution.
  • stabreg - Regression parameter estimates of a stable distribution.

The Tests and goodness-of-fit functions group includes only two functions:

  • edftests - Empirical distribution function (edf) goodness-of-fit statistics (Kolmogorov and Anderson-Darling).
  • ferrors - Compute price forecast errors (including MDE, MeDE and DRMSE).

The Markov regime switching (MRS) models group is composed of three complex functions for estimation, simulation and result visualization:

  • mrs_est - Estimate 2-state (Markov) regime switching model.
  • mrs_plot - Plot calibration results for a 2-state MRS model.
  • mrs_sim - Simulate trajectories of a 2-state (Markov) regime switching model.

The GUI functions group includes two utility routines:

  • onlyone - Sets only one active radiobutton.
  • readf - Read files and folders from the current directory.

The Demos group is probably the largest and most complex. It includes four main routines for illustrating the major Case Studies in the book, see Sections Distributions, Load forecasting, Price forecasting, Regime switching models for details. Each GUI routine makes use of many auxiliary functions:

  • mfe_distrib - GUI illustrating Case Study 2.6.3.
    Auxiliary routine(s): mfe_distrib_aux
  • mfe_loadf - GUI illustrating Case Studies 3.4.7 and 3.4.10.
    Auxiliary routine(s): mfe_loadf_block, mfe_loadf_crit, mfe_loadf_errtables, mfe_loadf_model_arma, mfe_loadf_plotacf, mfe_loadf_plotcomp, mfe_loadf_plotperiodogram, mfe_loadf_predA, mfe_loadf_predB
  • mfe_pricef - GUI illustrating Case Studies 4.3.7, 4.3.8 and 4.3.11.
    Auxiliary routine(s): mfe_pricef_errtables, mfe_pricef_plots, mfe_pricef_win_arx, mfe_pricef_win_arxg, mfe_pricef_win_tarx, startarx, startgarch, starttarx
  • mfe_regime - Script illustrating Case Study 4.4.5.

Finally, the Data files group includes four ASCII format datasets used in the Case Studies:

  • CA_daily.dat - daily electricity prices, loads, forecasted loads and max air temperatures from California (1998.04.01-2002.12.31).
  • CA_hourly.dat - hourly electricity prices, loads and forecasted loads from California (1998.04.01-2001.01.31).
  • EEX_daily_des.dat - deseasonalized (with respect to the weekly cycle) EEX daily price (Phelix Base index, 2002.01.01-2004.12.31).
  • NP_daily_des.dat - logarithm of the deseasonalized with respect to the weekly and annual cycles average daily Nord Pool system price (1997.01.03-2000.04.27).

Back to Contents


4. Distributions (Case Study 2.6.3)

To open GUI illustrating Case Study 2.6.3 type in the following command:

>> mfe_distrib

In return MATLAB will display the main window. After selecting the data file (use daily data; note, that the second column is used only) and pressing the 'Open/Load' button all objects will become active.

To preprocess the data by removing the 7-day cycle with the moving average technique (remst.m) tick the 'remove 7-day cycle' checkbox. To transform it further by taking first differences ('diff') or log-returns ('diff') select the appropriate radiobutton.

Next, choose which of the four distributions (Gaussian, hyperbolic, NIG, alpha-stable) will be fitted to the data and whether EDF test statistics (Kolmogorov and Anderson-Darling) will be computed. Finally, press the 'Run' button. The results will be displayed in the command window and visualized in 3-4 figures:


EEX deseasonalized price differences


Left tail of the distribution

Back to Contents


5. Load forecasting (Case Studies 3.4.7 and 3.4.10)

To open GUI illustrating Case Studies 3.4.7 and 3.4.10 type in the following command:

>> mfe_loadf

In return MATLAB will display the main window. It is divided into three thematic boxes. Controls of the top box allow for visualizing Model A (differencing - smoothing; for definition of N see formula (2.10) in [1]) and Model B (rolling volatility + moving average) seasonal decomposition approaches by plotting:

  • seasonal (deterministic) and stochastic components
  • periodograms of the stochastic components
  • ACF and PACF of the stochastic components

Controls of the middle box allow for finding the best - in terms of the AICC or BIC criterion - ARMA(p,q) model, with p = 1, ..., Max. AR order, q = 1, ..., Max. MA order, for the stochastic component of Model A or Model B. The results are displayed in the command window in two tables:

  • Best ARMA(p,q) model for a given AR order p = 1, ..., Max. AR order
  • Best 5 models

Finally, the bottom box allows for performing day-ahead forecasts of the 2001-2002 California system-wide load. 'AR', 'MA' and 'X order' edit boxes define orders p, q and r, respectively, of an ARMAX(p,q,r) model. Values of the exogenous variable (max air temperature) for days T-'X shift', T-'X shift'-1, ... are used to forecast load for day T. Afterwards a summary of the forecasting performance is displayed in the command window (MATLAB's format is used, i.e. ARMAX[p r q 'X shift']):

Back to Contents


6. Price forecasting (Case Studies 4.3.7, 4.3.8 and 4.3.11)

To open GUI illustrating Case Studies 4.3.7, 4.3.8 and 4.3.11 type in the following command:

>> mfe_pricef

In return MATLAB will display the main window. It has four buttons for selecting the forecasting model:

  • AR(X),
  • TAR(X),
  • AR(X)-GARCH,
  • AR(X) + GARCH,

and three edit boxes for defining the calibration period and the number of days to be forecasted.

If the 'ARX model' button is pressed the 'ARX model' window appears. The top box allows for:
  • specifying the order of the AR part,
  • selecting whether the exogenous variable (CAISO load forecast) is used ('-X' checkbox),
  • selecting the spike preprocessing method ('None', 'Similar day', 'Limit' or 'Damped', see [1], Section 4.3.8) and the preprocessing threshold coefficient (a spike is defined as any value exceeding 'mean(data)' + coeff * std(data)).

The lower box concerns interval forecasts and allows for:

  • specifying confidence levels (a maximum of three),
  • selecting whether 'gaussian' and/or 'non-parametric' intervals are computed (see [1], Section 4.3.12)
  • specifying whether interval (and point) forecasts are price capped
If the 'TARX model' button is pressed the 'TARX model' window appears. The top box allows for:
  • specifying the order of the AR part,
  • selecting whether the exogenous variable (CAISO load forecast) is used ('-X' checkbox),
  • selecting the threshold variable (price or load).

The lower box has the same functions as in the 'ARX model' window.

If the 'ARX-GARCH model' button is pressed the 'ARX-GARCH model' window appears. The top box allows for:
  • specifying the order of the AR part,
  • selecting whether the exogenous variable (CAISO load forecast) is used ('-X' checkbox),
  • specifying the orders of the GARCH(P,Q) part.

The lower box has the same functions as in the 'ARX model' window.

In the ARX-GARCH model parameters of the ARX and GARCH parts are estimated simultanously. Unfortunately, MATLAB's garchfit.m does not do a very good job of it - the price forecasts are very noisy. In [1] SAS was used to obtain the forecasts. The forecasts were better than MATLAB's, but still worse than those of the ARX model.

The 'ARX + GARCH model' button opens an almost identical window (except for the window name). In this model first the ARX part is estimated, then the residuals are modeled by a GARCH process. The resulting forecasts are better than for the ARX-GARCH model, but still slightly worse than those of the ARX model.

During the estimation/forecasting process, current status (forecast n out of N forecasts) is displayed in the command window. Afterwards a summary of the forecasting performance is displayed in two tables - one with daily errors (MDE, MeDE) and one with weekly errors (MWE, MeWE, WRMSE) for full weeks in the forecasting period:

Additionally the results are visualized in two figures - one with point and interval (here: 90% and 95% two-sided non-parametric) forecasts and one with deviations of point and interval forecasts from the actual market price:

Finally, the forecasting results are saved in an ASCII file finalresult.txt using the following format (a maximum of 29 columns in total):

  • date,
  • hour,
  • actual price,
  • price forecast,
  • capped price forecast,
  • [{lcl_gauss ucl_gauss, lcl_npar ucl_npar} price_capped{lcl_gauss ucl_gauss, lcl_npar ucl_npar}] as many times as there are confidence levels (but no more than 3),

where lcl stands for 'lower confidence level' and ucl for 'upper confidence level'.

Back to Contents


7. Regime switching models (Case Study 4.4.5)

To run the script illustrating Case Study 4.4.5 type in the following command:

>> mfe_regime

It loads Nord Pool deseasonalized log-prices (NP_daily_des.dat file) and estimates 2-state Markov regime switching (MRS) models with the base regime driven by a mean reverting process (i.e. AR(1)) and the spike regime modeled by one of the following:

  • Gaussian random variable,
  • lognormal random variable,
  • Pareto random variable,
  • mean reverting process (like for the base regime).

The obtained estimates are displayed in the command window:

and visualized, here for Pareto (left) and mean-reverting (right) models:

In the bottom panels the probability Pr(spike) that an observation comes from the spike regime is plotted. In the top panels the time series is displayed with observations identified as belonging to the spike regime (i.e. having Pr(spike)>0.5) plotted in red.

The toolbox also allows for calibration of a MRS model with Weibull distributed spikes. To see the results change 'G' to 'W' in [Ksi_tT, Param, P, Ksi_t1t_10, LogL] = mrs_est(Data, 'G'). However, with the Nord Pool dataset the estimates do not converge. For simulation of MRS models see mrs_sim.m.

Back to Contents


MFE webpage

Homepage


Last modified on 2009-10-06.