Mathematical Bureau
Mathematical Modelling In Electricity Market

## Introduction

Recently, I made an overview of The short term electricity price forecast model on the most similar pattern which I developed during 2008-2013. You're welcome to read it. If you need an example of the simplest version to try out for your short-term forecast problem, please do so.

This example aims to provide the model script which you can easily read, learn, understand, and further develop for the actual forecast problem you're hoping to solve.

The example file EMMSP_example_ENG.zip contains:

1. priceDataEurPriceZone.mat — dataset of the day-ahead market prices for the European price zone,
2. EMMSP_example_ENG.m — MATLAB script.

Once you download the archive, you should put both files in the same folder and run MATLAB. The example is done in MATLAB 2015b version; it will probably work on earlier versions as well.

Let's go through the script step by step.

## Initialisation

clear; clc;

The line allows us to clean the workspace and the command line.

load priceDataEurPriceZone priceDataEurPriceZone;

This line uploads the price dataset to the workspace. If you put the priceDataEurPriceZone.mat file in the same folder as the script file, there won't be any trouble. If not, then you should set up the folder for the dataset in command load().

ForecastMomentT = datetime('03.09.2012 23:00',...
'InputFormat','dd.MM.yyyy HH:mm');
P = 24;

As an example, we want to make a forecast at the time moment 03.09.2012 23:00. This means that the first forecast value refers to 04.09.2012 00:00. Forecast horizon value P defines the number of forecast values we'd like to get.

You may set any forecast moment and change the forecast horizon if you want to test the model. The dataset contains values from 01.09.2006 00:00 to 13.11.2012 00:00. Make sure you're stating the problem within this range.

% Forecast model parameters

M = 144;                % Model on the most similar pattern parameter
loopStep = 24;          % Additional variable for looping

Forecast model on the most similar pattern parameter M = 144. The value is obtained through the model calibration procedure. I do not discuss this procedure in this example.

## Model implementation

Descriptions for each step you may find here: The short term electricity price forecast model on the most similar pattern

### Step 1. Define the latest available actual pattern.

T = find(priceDataEurPriceZone.timestep == ForecastMomentT);

if ~isempty(T)

patternLatestAvailable.timestep = ...
priceDataEurPriceZone.timestep(T-M+1:T);
patternLatestAvailable.value = ...
priceDataEurPriceZone.value(T-M+1:T);

checkValues = patternLatestAvailable.value == 0;
if sum(double(checkValues)) == M
fprintf('Error: patternLatestAvailable contains only zeros',...
', check the data!');
return
end

else

fprintf('Error: forecast moment %s is absent',...
'in the time series, please check!', datestr(T));
return

end

This part of the script finds the piece of time series which precedes forecast moment — the latest available actual pattern. The time series piece is called a pattern when it comes to the forecast model on the most similar pattern. Also, I added a few error handlings.

### Step 2. Find the most similar pattern.

k = T - M - P;

similarityMeasure.timeDelay = [];
similarityMeasure.value = [];

% Looping through timeseries history
while k > M

patternTemp.value = priceDataEurPriceZone.value(k+1 : k+M);
patternTemp.timestep = priceDataEurPriceZone.timestep(k+1 : k+M);
similarityMeasure.timeDelay(end+1) = T - k;

checkValues = patternTemp.value == 0;
if sum(double(checkValues)) == M
similarityMeasure.value(end+1) = 0;
else
similarityMeasure.value(end+1) = ...
abs(corr(patternTemp.value,...
patternLatestAvailable.value,...
'type', 'Pearson'));
end

k = k - loopStep;       % looping not through hours, but through days

end

This while-looping is the most time-consuming part of the script. We're looping through all available patterns in the actual data and picking up the most similar one in the following lines.

maximumSimilarity = max(similarityMeasure.value);
msp = similarityMeasure.timeDelay(...
similarityMeasure.value == maximumSimilarity);
fprintf('maximumSimilarity = %1.3f',...
' with timeDelay = %d\n', maximumSimilarity, msp);

patternMostSimilar.timestep = ...
priceDataEurPriceZone.timestep( T-msp+1 : T-msp+M );
patternMostSimilar.value = ...
priceDataEurPriceZone.value( T-msp+1 : T-msp+M );

### Step 3. Calculate the linear coefficients.

A = regress(patternLatestAvailable.value,...
[patternMostSimilar.value ones(size(patternMostSimilar.value))]);

Linear coefficients are obtained with the least mean square method — just one line in MATLAB.

Step 4. Define the base pattern.

patternBase.timestep = ...
priceDataEurPriceZone.timestep( T-msp+M+1 : T-msp+M+P );
patternBase.value = ...
priceDataEurPriceZone.value( T-msp+M+1 : T-msp+M+P );

The base pattern is the one to be used to get the forecast values. For further explanation and illustrations of the pattern positions, please, check the blog on The short term electricity price forecast model on the most similar pattern.

### Step 5. Calculate forecast values.

To get the forecast values we merely need to apply linear coefficients to the base pattern.

Finally, we should estimate the error and plot the results.

patternActuals.timestep = ...
priceDataEurPriceZone.timestep( T+1 : T+P );
patternActuals.value = ...
priceDataEurPriceZone.value( T+1 : T+P );

MAE = mean(abs(patternForecast.value - patternActuals.value));
MAPE = mean(abs(patternForecast.value - patternActuals.value)...
./ patternActuals.value);

Firstly, I estimate forecast error by the two values: MAE and MAPE. MAE refers to the mean absolute error and MAPE applies to the mean absolute percentage error. These two values give us the first impression of the forecast quality. Furthermore, you can assess the number of different indexes and error dimensions for a better understanding of the model behaviour.

fprintf('Foreast error: MAE = %3.2f',...
' Rub/MWh, MAPE = %2.2f%%\n', MAE, MAPE*100);

plot(patternActuals.timestep,...
[patternActuals.value patternForecast.value])
title(sprintf('Day ahead price forecast',...
'forecast moment %s, MAE = %3.2f Rub/MWh, MAPE = %2.2f%%',...
datestr(T), MAE, MAPE*100))

At the very end of the script, we're printing the forecast error values and making the simplest version of the plot.

I repeat, this example can be uploaded in the file EMMSP_example_ENG.zip and run on your computer. A description of the time series forecast model on the most similar pattern is given in the blog The short term electricity price forecast model on the most similar pattern.