Estimation of Forward Curves

A First Expedition into the World of Electricity Prices

Preparation

<< MathExcel`    
<< Graphics`      
<< Statistics`   
<< TimeSeries`   
<< Finance`         
$DefaultFont  =  {"Arial", 14};   
Off[General::"spell1"]    
SetOptions[ListPlot, PlotJoined -> True];

The prices for electricity we want to analyze are provided from the Leipzig Power Exchange (LPX) as an Excel Sheet. The sheet can be downloaded from http://www.lpx.de/info_center/downloads/dl_general/LPX_Results.xls . To get started the file is opened with Excel and Mathematica is linked to the Excel sheet.

xl =  XLOpenLink[OpenMode -> omAttach]

A first look at the data

The first column within the Excel sheet contains the dates. These dates will come over the link as OLEDate(s). To get the typical Mathematica representation of dates these OLEDates are transformed to lists. (If we find time we should automate the selection of the range. At the moment the range is hard coded.)
In the Excel sheet the first data row is the most recent one. For time series analysis it is more convenient to have the most recent observation as the last element. The reordering of the data can be achieved with the Mathematica command Reverse. Since the first day in the historical data set starts with prices of zero for some unknown reasons, we just throw away this observation with the command Rest.

dates  =  Flatten[XLGetRange["A2:A685"]]  /.  
              OLEDate -> List  //  Reverse  //  
          Rest;
dates  //  Shallow

In the second column of the Excel sheet are the prices for the base load.

base  =  XLGetRange["B2:B685"] //  Flatten  //  
            Reverse  //  Rest;
base  //  Shallow

To get a first impression of the data we plot the data, doing some adjustments to the regular plot options and providing non-standard ticks.

ticks  =  
        Map[{#, StringForm["`3`/`2`/`1`", Sequence @@ dates[[#]]]} &, 
          Range[1, Length[dates], 100]];
plbase  = 
        ListPlot[base, PlotJoined -> True, AspectRatio ->  0.2, 
          Frame -> True, 
          FrameTicks -> {ticks, Automatic, None, None}, 
          PlotStyle -> Green, FrameLabel ->  "\:f280/MWh", 
          Axes -> None, PlotRange -> All];

Ther are large jumps observable in December 2001. At the moment we are not interested in these jumps therefore we use another scale in the plotting.

plbase  = 
        ListPlot[base, PlotJoined -> True, AspectRatio ->  0.2, 
          Frame -> True, 
          FrameTicks -> {ticks, Automatic, None, None}, 
          PlotStyle -> Green, FrameLabel ->  "\:f280/MWh", 
          Axes -> None];

There is no up- or downward trend recognizable in the data. Remarkable are the repeating patterns in the data. Observable are also jumps with different magnitude, whereby the level comes back very fast after a jump.

The third column within the Excel sheet has the prices for the peak load.

peak  =  XLGetRange["C2:C685"] //  Flatten  //  
            Reverse  //  Rest;
plpeak  = 
        ListPlot[peak, PlotJoined -> True, AspectRatio ->  0.2, 
          Frame -> True, 
          FrameTicks -> {ticks, Automatic, None, None}, 
          PlotStyle -> Red, FrameLabel ->  "\:f280/MWh", 
          Axes -> None];

The plot for the peak load looks at first sight very similar to the plot of the base load. The similarity will be analyzed later.

In the Excel sheets are also the prices for every hour.

hour  =  XLGetRange["L2:AI685"]  //  Reverse  // 
          Rest;
ListPlot3D[hour, BoxRatios ->  {2, 3, 1}, 
        Ticks -> {Automatic, ticks, Automatic}, 
        ViewPoint -> {2.848,  1.291,  1.294}, 
        PlotRange -> {0, 100}];

The warning from Mathematica tells us that some prices are missing. We will ignore this message for the moment. It looks like there are a lot of  jumps in the hourly prices. Also it looks like there is a typical pattern over a day. To further explore this pattern over a day the mean for each hour is calculated and plotted.

Mean[hour] 
BarChart[% /.  Null ->  0, 
        Frame -> {True, True, False, False}, 
        FrameLabel ->  {"hour", ""}, AspectRatio ->  0.2, 
        Axes -> None];

Comparing base and peak load prices

To further explore the seemingly similarity of the peak and base load prices both time series are plotted in one graph together with the differnces between them. As expected the prices of the peak load are higher then the base load in most cases.

DisplayTogether[Show[plbase, plpeak, PlotRange -> {0, 100}], 
        BarChart[peak - base, AspectRatio ->  0.2, 
          Ticks -> {ticks, Automatic}]];

But there are some exceptions to this rule

Count[peak - base, x_ /; x < 0]

The correlation which can be sensed from a look at the data is confirmd through a regressionen analysis. For the regression analysis we look at the one day changes of the prices.

regdata  = {Difference[peak, 1], Difference[base, 1]}  //  
          Transpose;
ListPlot[regdata, PlotJoined -> False, Frame -> True,  
        PlotRange ->  Automatic];
Regress[regdata, {1, y}, y]

Interesting is always the look at the residuen of the regression.

res  =  
        Regress[regdata, {1, y}, y, 
          RegressionReport -> { FitResiduals, 
              EstimatedVariance}];
ndist  =  
        NormalDistribution[0, (Sqrt[EstimatedVariance])]  /.  
          res;
resdata  = FitResiduals   /.  res;
DisplayTogether[Histogram[ resdata , HistogramScale -> 1] , 
        Plot[PDF[ndist, y], {y, -6, 6}]];

A mean to explore through visualization the possible deviation of the observed distribution to the normal distribution is the QQ-plot. The QQ plot is not available in Mathematica but it can be constructed with a few lines of Mathematica code.

\(QQPlot[data_, dist_, opts___?OptionQ] := 
      Module[{orderstats = Sort[data], n = Length[data], refQuantile, yTicks, 
          nline}, 
refQuantile\  = \ 
          Function[x, Quantile[dist, x], 
            Listable]; 
yTicks\  = \ 
          Map[{refQuantile[#/100. ], #} &, {1, 5, 10, 25, 50, 75, 90, 95, 
              99}]; 
nline\  = \ 
          Map[{refQuantile[#], refQuantile[#]} &, {0.01, 
              0.99}]; 
Show@
Graphics[{{SteelBlue, 
MapThread[
                  Circle[{#1, #2}, Offset[{3, 3}]] &, {orderstats, 
                    refQuantile[Range[n]/\((n + 1. )\)]}]}}, opts, 
            Frame -> 
              True, 
FrameTicks -> {Automatic, yTicks, 
                None, None}, FrameLabel -> \ "Quantiles", 
            Prolog -> {{GrayLevel[0], \ Line[nline]}}, 
            AspectRatio -> Automatic]
]\)

Applying the QQ plot to the residuals shows us that the distribution of the residual is not a normal distribution. If it would be, the blue circles would be closer to the straight black line.

QQPlot[resdata, ndist, 
        FrameLabel ->  "Quantiles"] ;

Looking for auto correlations

We are also interested in possible auto-correlations in the data. To visualize them we can build an autocorrelation plot, which is constructed around the function CorrelationFunction which comes with the package TimeSeries.

AutoCorrelationPlot[data_, k_, opts___?OptionQ]  :=  
      Module[{n = Length[data]}, 
        Show@
Graphics[{{SteelBlue, 
                Sequence @@ 
                  MapThread[{Disk[{#1, #2}, Offset[{3, 3}]], 
                        Line[{{#1, 0}, {#1, #2}}]} &, {Range[k + 1] - 1, 
                      CorrelationFunction[data, k]}]}, {Dashing[{0.01}], 
                GrayLevel[
                  0.5], 
Line[{{0, 2/Sqrt[n]}, {k, 
                      1.96/Sqrt[n]}}], 
Line[{{0, -2/
                        Sqrt[n]}, {k, -1.96/Sqrt[n]}}]}}, opts, 
            Axes -> {True, False}, Frame -> True, 
            FrameTicks -> {Automatic, Automatic, None, None}, 
            PlotRange -> All, 
            FrameLabel -> {"Lag",  ""}]]

The doted gray lines are the boundary at which autocorrelations are significant at a 95% level.

Day-of-week effects

The newly constructed plot can be used for analyzing the data. As a first analysis the one day differences in the data are explored. The function Difference comes from the TimeSeries package.

AutoCorrelationPlot[Difference[peak, 1], 30, 
        AspectRatio -> 0.3];

There seems to be significant autocorrelation in respect to the previous day, the day before the previous day and then on the seventh, fourteenth, ... (7^n) th day. This indicates a day-of-the-week effect in the data.

Seasonal effects

If we only look on weekly sequences of the data a seasonal pattern can be observed.

DayOfWeek[dates[[3]]]

AutoCorrelationPlot[peak[[Range[3, Length[peak], 7]]], 52, 
        AspectRatio -> 0.3];

When looking at all the sundays in the data there is a significant negative autocorrelation around the 26th week, indicating a half year effect. It is interessting that these autocorrelation are not so clearly visible for other days of the week (yet?)

DayOfWeek[dates[[6]]]

AutoCorrelationPlot[peak[[Range[6, Length[peak], 7]]], 52, 
        AspectRatio -> 0.3];

A first try to construct a forward curve

We have the feeling that parts of the price changes are not random, i.e we think we are able to predict in part the future evolvement of the prices. To achieve this we try to identify the non-random parts of the data and build a function for them. We call this function a forward curve. Our assumption is that the future prices will randomly realize around this forward curve.

Adjusting for weekly pattern

Doing an autocorrelation analysis and looking at the data seemed to reveal an day-of-the-week effect. To further explore this day-of-the-week effect, the data is rearranged in weekly chunks, then the mean for each week is calculated and the difference of each day of the week to the mean of the week is calculated.

data  =  
        Map[Transpose[{Range[7], # - Mean[#]}] &, 
            Partition[Drop[peak, 3], 7]]  //  Flatten[#, 1] &;

This results in the following picture. The x-axes describe the week days starting with Monday and ending with Sunday. As expected are the prices on Sunday below the average of the week.

plt  = 
        ListPlot[data, PlotJoined -> False, 
          PlotStyle -> {Green, PointSize[0.02]}, Frame -> True, 
          PlotRange -> Automatic];

To fit a curve into this data we use our usual first template for non-linear estimation a + 
        b (e^(s x)).

wfkt[x_]  = 
      NonlinearFit[data, a + b Exp[s x], x, {a, b, s}]

In this case the function seem to able to follow the data.

fktplt  =  
        Plot[wfkt[x], {x, 1, 7}, 
          PlotStyle ->  {Thickness[0.01], Yellow}, 
          DisplayFunction -> Identity];
Show[plt, fktplt];

We then use the function to "correct" the observed prices. The RotateLeft is necessary because the data does not start with a Monday.

peak1  =  
        peak  - 
          Map[wfkt[#] &, 
            RotateLeft[Mod[Range[Length[peak]], 7] + 1, 3]];

The autocorrelation plot for the "corrected" data reveals that the weekly autocorrelation has disappeared.

AutoCorrelationPlot[Difference[peak1, 1], 30, 
        AspectRatio -> 0.3];

Controlling for holidays

From the observation, that Sundays have special prices one might also conclude special prices on holidays. To explore this idea first we use the HolidayQ function which comes with our package Finance, which was loaded in the beginning. We do not have a holiday calendar for Leipzig available so we take Frankfurt instead.

hol  = 
        Position[dates, x_ /;  HolidayQ["FRA", x]]  //  
          Flatten;
dates[[hol]]

Marking the holidays in the time series with a green point convinces us that our conclusion about holidays is not false:

DisplayTogether[
        ListPlot[peak1, PlotJoined -> True, AspectRatio ->  0.2, 
          Frame -> True, PlotRange ->  {0, 70}, 
          FrameTicks -> {ticks, Automatic, None, None}, 
          PlotStyle -> Blue, FrameLabel ->  "\:f280/MWh", 
          Axes -> None], 
        ListPlot[Transpose[{hol, peak1[[hol]]}], 
          PlotJoined -> False, 
          PlotStyle -> {PointSize[0.015], Green}]];

We calculate the relative mean price decrease due to a holiday.

prev  = 
      Mean[((peak1[[hol - 1]] - peak1[[hol]])/peak1)[[hol - 1]]]

and the relative price increase after a holiday.

next  = 
      Mean[((peak1[[hol + 1]] - peak1[[hol]])/peak1)[[hol + 1]]]

The simplest approach controlling for holidays is to get rid of them. Therefore we throw them away

peak2  = Delete[peak1, Map[List, hol]];
dates2  =  Delete[dates, Map[List, hol]];

Estimating seasonal effects

In the autocorrelation analysis were seasonal effects visible at least when looking on the Sunday's prices. Other studies of electricity prices from other countries also regulary find seasonal pattern. To estimate seasonal patterns the dates are expressed as fractional years.

t  = YearsBetween[{2000, 1, 1}, dates2];
t  //  Shallow

and the prices are calculated relative to the mean.

m  =  Mean[peak2];
data  =  
        MapThread[{#1, #2 - m} &, {t, peak2}];
data  //  Shallow

For estimating a functionf for the seasonal effects a non-linear curve is fitted on this data. The template for the function was taken from Pilipovic (1997).

yfkt[y_]  = 
      NonlinearFit[data ,  a y + b Cos[c π y], 
        y, {a, b, c}]

To compare visually the estimated curve with the data the following graph is constructed

ticks2  =  
        Map[{#, StringForm["`3`/`2`/`1`", 
                Sequence @@ dates2[[#]]]} &, 
          Range[1, Length[dates2], 75]];
DisplayTogether[
        ListPlot[peak2, AspectRatio ->  0.2, Frame -> True, 
          FrameTicks -> {ticks2, Automatic, None, None}, 
          PlotStyle -> Red, FrameLabel ->  "\:f280/MWh", 
          Axes -> None, 
          PlotRange -> Automatic], 
ListPlot[
          m  + Map[yfkt, t], 
          PlotStyle -> {Thickness[0.005], Blue}]];

Fur further analysis (not shown in this notebook) the estimated function for seasonality is used to further adjust the observed prices.

peak3  =  peak2  - Map[yfkt, t];

Constructing the forward curve

Using the results of the previous analysis a forward curve can be constructed for the valuation date. It is assumed that the valuation date is the next to the last observed date. The forward curve should start at the last observed price, assuming that the last price together with the above estimated corrections is the best estimation for the future prices.

valdate  = DaysPlus[Last[dates2], 1]
DayOfWeek[%]
(Subscript[S, 0]) =  Last[peak3]

The dates for the coming 365 days are constructed and the day of the week and the year fraction for each day are calculated.

dates3  = 
        Map[DaysPlus[Last[dates2], #] &, 
          Range[1, 365]];
week  = Mod[Range[365], 7] + 1;
year  =  YearsBetween[{2000, 1, 1}, dates3];

Together with the above estimated curves a forward curve can be constructed.

forwards  =  
        Transpose[{dates3, year, 
            (Subscript[S, 0]) + Map[wfkt, week] + Map[yfkt, year]}];

To visualize the estimated forward curve it is added to the plot of historical prices.

n  =  Length[peak];
ticks3  =  
        Map[{n + #, 
              StringForm["`3`/`2`/`1`", Sequence @@ dates3[[#]]]} &, 
          Range[1, Length[dates3], 100]];
plforwd  = 
        ListPlot[Transpose[{Range[n + 1, n + 365], forwards[[All, 3]]}], 
          PlotStyle -> Green, 
Frame -> True, 
          FrameTicks -> {ticks3, Automatic, None, None}, 
          DisplayFunction -> Identity];
Show[plpeak, %, PlotRange -> Automatic, 
        FrameTicks -> {Join[ticks, ticks3], Automatic, None, 
            None}];

Summary

This example for estimating forward curves are a first expedition into the world of electricity prices. The analysis was done in preparation for a project. Doing this exercise helped us to understand more sophisticated approaches for estimating forward curves. Several parts are missing in this notebook like the estimation of volatilities, adjusting the forward curve to trade prices, controlling for jumps etc. Eventhough we left out a lot of important topics  we hope that this example is usefull as a starting point for others.

References

Cleveland, William S., (1993), Visualizing Data, New Jersey 1993.

Clewlow, Les/Strickland, Chris, (2000), Energy Derivatives: Pricing and Risk Management, Lacima Publications.

Pilipovic, Dragana, (1997), Energy Risk: Valuing and Managing Energy Derivatives, McGraw-Hill.